#### Abstract

With the development of society, public resources for healthcare are increasingly inadequate to meet the demands for the services. Therefore, it is extremely important for policymakers to provide citizens with the most effective healthcare services within the limited available resources. In order to achieve positive effect rescue operations in the Emergency Medical Services (EMS) system, the problems including where to locate the ambulance facilities and how many ambulance vehicles should be allocated to the stations have become the focus of attention. In this paper, we study the problem based on the demand for EMS in Songjiang District, Shanghai, China, followed by the joint planning of Emergency Medical Services management, which typically consists of ambulance facility locations planning and patient’s assignment to hospitals. We proposed a modified Double Standard Model (DSM) to maximize the demand points covered at least two times within the minimum coverage criteria. The problem is solved by integer linear programming technique with the CPLEX software and we make a comparison between the solutions and the locations which exist in the emergency system used by the Songjiang emergency center. Our results show that the demand coverage rate and response time can be efficiently improved through relocating the current facilities without additional vehicle resources.

#### 1. Introduction

Because the Emergency Medical Services (EMS) provide the fundamental assistance and transportation for patients in need of urgent medical care in saving lives, numerous studies have been carried out in order to improve the quality of the EMS. However, under the pressure of the financial budget burden to meet the standard levels of performance criteria, EMS decision makers have to deal with the difficult tasks: where to locate the ambulance facilities, how many emergency vehicles can be allocated to the ambulance stations, and so forth. These kinds of problems become more and more complicated with the arrival of the aging society and the high population mobility, especially serious in the cities such as Shanghai and Beijing in China. According to the statistical data from Shanghai Municipal Commission of Health and Family Planning (SMCHF) in the year of 2014, the emergency calls in Shanghai climbed up to more than over 642 thousand times.

Among the EMS literature, the problems of allocating ambulance vehicle resources and the location decisions of ambulance station have attracted special attention during the past decade of research (cf. Brotcorne et al. [1]; Li et al. [2]; Aringhieri et al. [3]). In order to improve the emergency efficiency for the service demand point, EMS departments manage to locate more ambulance vehicles or to add more ambulance stations. However, it is not reasonable to use a large number of vehicles or stations to cover the emergency demands. The new ambulance station construction cost is high and the ambulance vehicles are expensive (Su et al. [4]). Each service vehicle is completely equipped with all necessary emergency facilities that medical team may need in their missions (Moeini et al. [5]). Hence, with a limited number of ambulance vehicles, we need to have a balance between the ambulance stations and the desired objective.

Though there is a substantial body of work on locating the emergency service vehicles and facilities planning, most existing models, solution approaches and studies in this research area concentrate on planning without considering the relationship between the ambulance facility locations and the hospitals destination. As one of the important aspects in the EMS, the selection of hospitals also plays a crucial role in the Emergency Medical Service for patients. Generally speaking, most studies considered that early detection, early response, and appropriate medical assistant are important factors in the EMS system. However, they are just the first stage in Emergency Medical Service. We should consider the distance between the demand point and the hospital, patients’ injury degree, and the ability and capacity of rescue team in hospitals as the second stage in the EMS system. Traditionally, most patients are sent to the nearest hospitals. As we all know, there is a hierarchy for the services provided by hospitals. In China, there are three different medical levels, where each level in a hospital has its own medical resources. On one hand, some badly wounded patients must be sent to the advanced hospitals; on the other hand, advanced hospitals have their own capacity limits, and not all the patients can be served by the advanced hospitals. For the sake of providing the first aid to patients efficiency in the emergency system, decision makers need to consider the demand characteristics, such as patients’ injury levels, and the existing hospitals’ service levels. To the best of our knowledge, only Yang et al. [6] have considered the two-stage emergency rescue in the EMS with the objective of the total social cost minimization, but they did not consider the ambulance vehicle busy fraction in their model.

Under the above consideration, we are interested in proposing a new mathematical model that investigates the joint planning of ambulance station location, ambulance vehicle allocation, and the assignments of the ambulance stations to demand points and demand points to the hospitals. The new model is based on the Double Standard Model, called DSM. The DSM has been introduced by Gendreau et al. [7]. It aims to maximize the portion of demands being doubly covered within a primary standard and also ensures a full coverage within a secondary response standard. Both service coverage standards are measured as the travel time in minutes required by an ambulance vehicle to travel from an ambulance station to a demand point. In the DSM, all the patients’ injury level is implicitly assumed to be uniform. Therefore, patients cannot deal with the selections of ambulance stations and the corresponding hospitals. Also, in the DSM, the double service coverage treats each demand point separately and uniformly, which may generate results of service capacity redundancy for some demand points and shortage for other points. The chance constraint approach establishes a performance based reliability level by simultaneously considering the EMS vehicles and competing demand points, which could avoid the above limitations of imbalanced EMS vehicle allocation from the multiple service coverage approach (Liu et al. [8]).

The remainder of the paper is organized as follows. In Section 2, a brief literature review is given. Section 3 states the studied problem, and a mathematical formulation is proposed in Section 4. The model is tested on real-world data sets collected through the EMS system of Songjiang District, Shanghai. The computational experiments are reported in Section 5. Finally, the last section summarizes the work and indicates further research directions.

#### 2. Related Work

Emergency Medical Service is a classical problem; it has attracted researchers’ attention in the last decade. There is a stream of literatures dedicated to the study of EMS, as it is vital and plays a preliminary function in providing high quality healthcare services in our life, which could give the first answer in distressing situations. We give a brief literature review which includes the most related references on ambulance location and the hospital selection.

##### 2.1. Ambulance Location (Deterministic Models)

Most early proposed ambulance location models were integer linear formulations. Since these models did not consider the probability that an ambulance might be busy at a given time, they were classified as deterministic (Rajagopalan and Saydam [9]). The earliest EMS models have been introduced in the 70s by the articles of Toregas et al. [10] and Church and ReVelle [11]. Toregas et al. [10] studied the Location Set Covering Problem (LSCP) which identifies the minimum number of facilities and their locations that cover all demand points within a certain distance. However, the LSCP model cannot cover all demands with limited resources in reality, and one important thing is that once an ambulance is dispatched, some demand points are no longer covered. Due to the fact that LSCP model treats all demand points identically, the solution may require more ambulances than actually needed or underestimate the number of ambulances needed for those locations with relatively heavy demand (Saydam and Aytuǧ [12]). To avoid such limits, Church and ReVelle [11] proposed the Maximal Covering Location Problem (MCLP), which locates a fixed number of facilities so as to maximize the amount of demand that is covered by at least one facility.

In each of the above models, a common problem is that coverage may become inadequate when some ambulance vehicles are busy. Hence, to compensate this shortcoming, most literature about the location problem of EMS have followed these two early studies by Toregas et al. (LSCP model) and Church and ReVelle (MCLP model). For example, Daskin and Stern [13] proposed a Hierarchical Objective Set Covering Problem model (HOSC), a hierarchical model for the LSCP, with the objective of minimizing the number of ambulance locations providing full coverage within a distance standard first and then maximizing the number of demand points under multiple coverage. As the HOSC may privilege the congestion of the ambulance, Hogan and ReVelle [14] proposed the Maximal Backup Coverage models BACOP 1 and BACOP 2. These models use two ambulances to cover the demands. Gendreau et al. [7] developed a model known as DSM in which all demands must be covered by ambulance located in a secondary coverage radius of minutes, and in addition, a certain proportion of the demand must be covered in a primary coverage radius of minutes. Liu et al. [8] proposed the service reliability for the demand points to make sure of double coverage.

##### 2.2. Ambulance Location (Probabilistic Models)

Deterministic models did not account for the probability of a particular ambulance being busy at a given time and overestimated the actual coverage provided [9]. Hence, to compensate for this shortcoming, probabilistic models were later developed. Probabilistic location models recognized that any given ambulance may be busy when it is called. Such uncertainty can be modeled via programming techniques or using a queuing framework (Rajagopalan et al. [15]). One of the first probabilistic models for ambulance location is the Maximum Expected Coverage Location Problem (MEXCLP), which appears in Daskin [16]. He removed the implicit assumption of deterministic models; that is, all units are available at all times, by assuming each ambulance has the same probability , called the busy fraction, of being unavailable to answer a call. He also assumed all ambulances are independent. Given a predetermined number of response units, MEXCLP maximizes the expected coverage subject to a response time standard.

##### 2.3. Hospital Selection

For the hospital selection problem, a widely used rule consists of selecting the closest hospital to the demand points (Su and Shih [17]; Aboueljinane et al. [18]). However, Ingolfsson et al. [19] mentioned that the proximity could not be the only influencing factor as in Edmonton EMS, as 50% of patients are not transported to the closest hospitals. Savas [20] described some other factors that should be considered, such as the available capacity of hospitals, the hospital having the appropriate specialists, particular equipment for the patient, and the patient’s preference choice. Repoussis et al. [21] proposed a mixed integer linear programming formulation combining ambulance dispatching, patient to hospital assignment, and treatment ordering problem in mass-casualty incidents. Also, some researchers used the empirical distribution as estimated from historical data to select the area of the destination hospital (Ingolfsson et al. [19]; Repede and Bernardo [22]).

With respect to the reviewed literature mentioned above, the requirements of locating ambulance facility in Emergency Medical Services are not onefold, and all the factors should be dealt with collectively, for practical purposes. Consequently, the main contribution of our study is that we propose a new mathematical model to address the problem of ambulance facility location jointly with hospital selection in the situation where there are patients with different types of injuries.

#### 3. Problem Statement

In this section, we first present the concerned problem in detail and then list the main assumptions.

##### 3.1. Problem Description

Similar to the events sequence discussed by Nogueira et al. [23] for the EMS calls, Figure 1 depicts the sequence of events associated with an EMS management. An ambulance system begins immediately after an accident, such as traffic, natural disaster, and rebellion, and finishes once the patient has been delivered to a hospital. Some disease’s golden period for rescue is fairly short, such as heart disease in which a patient needs to receive help typically less than five minutes. The EMS department managers usually consider how to locate the ambulances in a coverage to access patients quickly. Most exciting literatures have considered the similar problems. Especially in 1973, the United States Emergency Medical Services Act set the “coverage” standards as follows: in urban areas, 95 percent of calls must be reached in 10 minutes, and in rural areas, calls should be reached in 30 minutes or less [24]. However, this is just the first phase in the EMS, as when an ambulance arrives at the demand point, only some simple treatments can be done, such as bandaging and CPR. This is because the ambulance can only load some basic medical devices and the patients’ diseases are various and complex. Therefore, how to transport a patient to a hospital facility quickly is also a curial decision considered in the EMS. Patients injury type is the first important factor considered in the transportation, and the hospital treatment capacity is another important factor which determines whether the patient can be transferred and accommodated. We call it as the second stage in the EMS when transferring the patient to the hospital.

In this paper, we propose and formulate a joint location planning and ambulance assignment model based on the classic DSM. The objective of this model is to maximize the total number of demands covered by the primary coverage twice. Let be the set of potential vehicle locations, and let be the set of hospitals. The set of demand locations is denoted by . The problem is defined on two graphs, and , where is the set of edges and is the set of edges . With each edge , the travel time from node to node is , and with edge (), the travel time from node to node is . Each demand location is associated with a corresponding demand value . The variable denotes the demand patient point can be served by ambulance station and be sent to hospital . All demand locations must be reached within minutes and at least have a proportion of of all the demand that must be covered by minutes. The largest number of available vehicles can be located at a station denoted by .

When an emergency call is received from a demand point , the EMS department should send an ambulance vehicle to service it. But the ambulance may be busy, serving a competing demand at that time. As such, one or more vehicles in the neighborhood of the demand point that could meet the service coverage standard should be targeted for the deployments to provide reliable service. In order to estimate the lowest number of ambulance vehicles needed to be deployed around demand point to ensure a certain level of service reliability, a modified busy fraction of point , denoted by , is introduced by [8] as follows: in which the given parameter is the average service time by ambulance vehicle, the fixed parameter is the frequency of requests for emergency calls by competing demands around demand point per day, denotes the set of ambulance station locations within the coverage of demand point , denotes the set of competing demand points around demand point , denotes the ambulance number in station , and is regarded as a given ratio. Assume the emergency calls are independent, and all demand points have the same local busy fraction, each following a binomial distribution. So the probability of having at least one vehicle available can be computed as follows:

With a predefined service reliability level, , for demand point , to meet the requirement of having at least one ambulance vehicle that can service it within the coverage , this constraint can be expressed as

The total number of ambulances in station can be expressed as follows:where is the smallest integer meeting the following expression:

##### 3.2. Assumptions for the Study

Some assumptions for this deterministic facility location and ambulance allocation problem are made as follows.(1)We consider all hospitals , classified as level I (the higher level of care) and level II (the lower level of care).(2)Each hospital has a fixed capacity of , treating patients.(3)To make full use of the hospital resources and give effective first aids to the patients, we also classified all the patients into two categories, as level I (the higher level of damage) and level II (the lower level of care).(4)Once the primary treatment on scene is completed, a destination hospital has to be selected. Patients in damage level I can only be sent to the hospital level I, and the damage level II can be sent to either hospital level I or hospital level II.(5)In our research, the traveling time from the demand point to the hospital should not exceed thirty minutes and the patient should be sent to the nearest hospital with appropriate service.

#### 4. Mathematical Model

In Section 4.1, we formally give the definitions of the problem parameters and the decision variables. In Section 4.2, we present a mathematical model for the problem.

##### 4.1. Notation

*Index Sets* : set of the demand points indexed by , . : set of the potential ambulance locations indexed by , . : set of the hospitals indexed by , . : set of ambulance station locations with coverage of the demand point with . : set of competing demand points around demand point , with any ambulance station and .

*Problem Parameters* : upper bound of the total ambulance vehicles. : upper bound of the ambulance in each station. : the type of the hospital : 1 means it has advanced medical level and 2 means otherwise. : the type of the patient in demand : 1 means patient is badly hurt and 2 otherwise. : the service reliability level. : the treatment capacity of hospital . : the density of the emergency demand at point . : the demand point ’s ambulance vehicle busy fraction. : the portion of demand points covered within the primary standard . : the smallest number of vehicles assigned around demand point at ambulance station , which satisfies . : average service time including the time from ambulance station to demand point, demand point to hospital plus the time from hospital to ambulance station, and the corresponding round trip of an ambulance. : time needed for ambulance from ambulation station to demand point . : time needed for ambulance transferring patient from demand point to hospital . : a binary parameter, equal to 1 if minutes (secondary standard), 0 otherwise. : a binary parameter, equal to 1 if minutes (primary standard), 0 otherwise. : a binary parameter, equal to 1 if minutes, 0 otherwise. : a binary parameter, equal to 1 if patients’ damage level matches hospital ’s level, that is, , 0 otherwise. : a sufficiently large integer.

*Decision Variables* : binary, equal to 1 if ambulance at station services demand point and sent the patient to hospital , 0 otherwise. : binary, equal to 1 if ambulance station is open, 0 otherwise. : binary, equal to 1 if demand point is covered one time by an ambulance vehicle within primary standard at reliability level , 0 otherwise. : binary, equal to 1 if demand point is covered two times by an ambulance vehicle within primary standard at reliability level , 0 otherwise. : integer, the number of vehicles located at station .

##### 4.2. The Formulation of the Problem

Objective (6) means to maximize the demands in primary coverage covered two times. Constraint (7) means that every demand must be serviced by one ambulance station and sent to a hospital. Constraint (8) means that every demand must be serviced by at least one ambulance vehicle at least. Constraint (9) ensures that a proportion of the demands must be covered in primary standard .Constraint (10) ensures that any demand point can only be serviced by an opened ambulance station. These constraints link variables and . Constraint (11) means that the number of ambulance vehicles in each ambulance station cannot exceed . Constraint (12) guarantees that if the ambulance station is constructed, there is at least one ambulance vehicle. This constraint links variables and . Constraint (13) denotes that the total number of ambulance vehicles should not exceed . In constraint (14), the left-hand side is the total number of the vehicles that covered site within primary standard , and it should not be less than the summation of single and double coverage indicators, which is the right-hand side of the inequality. Constraint (15) ensures that demand site can be doubly covered only if it is already singly covered. This constraint links variables and . Constraint (16) is a service reliability constraint, which ensures that point must be covered by ambulances with service reliability level . Constraint (17) ensures that the total number of patients transferred to a hospital should not exceed its own capacity. Constraints (18) ensures that the travel time or distance between demand site and the hospital should not exceed minutes. Constraint (19) guarantees that the patient injury type must be matched with the hospital level.Constraint (20) defines the domains of the decision variables.

#### 5. Case Study and Computational Results

In this section, we apply our model to the ambulance assignment optimization based on data obtained from Shanghai Trauma Emergency Center.

##### 5.1. Database Analysis

Our case study corresponds to the EMS system in Songjiang District, Shanghai, China. Songjiang has a size of about 605 km^{2} and a resident population of approximately 1.75 million (2015). The urban extension is about 23.8 km in the north-south direction, 24.6 km in the east-west direction, and 39.3 km in the southwest-northeast direction. It currently has 17 administrative regions and the average population density is about 2893 inhabitants per square kilometer.

###### 5.1.1. Demand Points

The large number of EMS call locations makes it difficult to handle the ambulance station location optimization problem. In order to simplify the complexity of the ambulance location problem, most previous studies considered the center point of the administrative regions as the demand points [6]. But for some cities, the boundary of administrative regions is irregular, and thus it is hard to find the central points. And even worse, for some small countries, the administrative regions’ boundary is fuzzy. Based on the above consideration, we organize the county into 2-kilometer 2-kilometer square nodes with Arcgis 10.2 software. This gives us a total of 187 demand nodes and we assumed that the ambulances can be posted in any one of the nodes except in the nodes that make up the boundary. This can be seen in Figure 2 which displays all the demand points. Each demand point has patients of one or both of badly injured type and generally injured type. The EMS call center of the city receives more than 25,000 calls in the year of 2013, but just about 20,000 calls require ambulance vehicles. The spatial distribution of the whole year in 2013 can be seen in Figure 3. And the number of two types of patients in each demand point is seen in Figure 4.

###### 5.1.2. Demand Scenarios

To reflect the realistic demands, EMS calls would be further classified based on their time and date. Days are classified into workdays (Monday through Friday) and off-days (weekends and holidays) and the call times are arranged into four periods of a single day: (i) dawn (00:00 AM–6:00 AM), (ii) morning (6:00 AM–12:00 PM), (iii) afternoon (12:00 PM–18:00 PM), and (iv) evening (18:00 PM–24:00 PM). The number of emergency calls of the patients in each of these 8 scenarios is shown in Table 1.

###### 5.1.3. Potential Ambulance Station Locations and the Hospital Locations

We selected a total of 17 potential ambulance location sites in the city of Songjiang including the hospitals and community health service centers. The current 8 ambulance stations are also in our potential location sites and we add the remaining potential stations in each administrative region mentioned in the previous context. All the potential ambulance stations’ information in detail can be seen in Figure 2 and Table 2. From the detailed emergency call records, we find that there are more than 50 hospitals that are the potential destination hospitals. But some hospitals received the amount of patients in single digit in the whole year. Finally, we select 29 hospitals as the destination hospitals, each of which has received more than ten patients according to the emergency call records. Table 3 shows the detailed coordinates data for these hospitals. All the geographic coordinates of these potential stations and hospitals are obtained from the Baidu API with JAVA program; then we convert them into WGS-84 coordinates through Arcgis 10.2 software.

##### 5.2. Results and Discussion

In our experiments, we suppose that each of the emergency calls must be covered in less than 12 minutes; that is to say, the secondary coverage is twelve minutes following the standards set from the 13th Five-Year Plan of Shanghai prehospital service of first aid by the Shanghai Municipal Commission of Health and Family Planning (SMCHF) in the year of 2016. At least of the calls must be covered in 8 minutes with the prespecified service reliability level of , which means and . The number of ambulances in the system is , which is the exact number of ambulances operated in 2013. Furthermore, the maximum number of ambulances per station is 4.

The shortest travel time from ambulance stations to demand points and from demand points to hospitals is calculated on the basis of Euclidean distance. The running average speed of an ambulance is 55 km/ph, which is consistent with the traffic condition for most of the Shanghai urban streets.

Since in the proposed model is determined by the ratio and service reliability level using (5), it is imperative to compute by using (1) prior to estimating . The service call frequency at each demand point uses data in the year of 2013 calls, and we get the day average call per demand point. The average service time for EMS vehicles is set to be 0.9 h, which is calculated through 19,447 EMS calls received in the year of 2013. With the above parameter values, the steps to compute ratio for a demand point are as follows [8]: (i) finding all the ambulance station sets which can service demand point , that is, ; (ii) determining all demand points, including , that are within the coverage of 5 minutes for each ambulance station in , then summing their competing demands , and ensuring no duplication; (iii) computing the ratio by (1). Finally, the smallest integer can be found by iterations (or binary search).

The solutions reported in the next section have been carried out on Intel Core™ i7-6700K CPU @ 4.00 GHz with 32 GB of RAM. The standard solver Gurobi (version 7.0.1) has been used to solve the integer linear programming model, coded with C++. Since the size of the model is not too large, solution time for each integer program run is less than ten minutes.

###### 5.2.1. Current Ambulance Station

We start from investigating the performance of the 8 stations and 29 ambulances currently utilized by the EMS system in Songjiang District. We attempted to solve our model by limiting potential bases to these 8 stations. The resulting optimization problem is infeasible, indicating that the current bases cannot satisfy the predefined EMS standards. This is the case even if the secondary coverage fraction is set to be 0, which indicates that the current bases cannot ensure the second coverage within minutes. In fact, the smallest value of for which we could obtain a feasible solution was 13 minutes. Even so, feasibility is attained only for a maximum value of = 0.3, which corresponds to 30% coverage within minutes. The conclusion is that the current bases are far from achieving both primary and secondary coverage requirements. Those two cases are listed in Table 4.

###### 5.2.2. Is There a Feasible Solution of 8 Ambulance Stations?

Having failed to meet the predetermined standards with the current 8 bases, we seek to determine whether any other set of 8 bases would satisfy those requirements. In order to handle this problem, we need to add an additional constraint on the maximum number of bases. To find out the reasons why the optimization based on current ambulance stations fails to meet the coverage standard, we used the same quantity of ambulance stations. After solving the above model, we found that it is possible to locate a different set of ambulances from all the potential ambulance stations which fulfill the full coverage within minutes, with the limitation. When we used the problem 2’s standard minutes, we can get the proportion with . Table 5 shows the trade-off between the proportion and the objective value. Figures 5 and 6 show the 8 optimal ambulance locations within all the demand points satisfying the service coverage in 13 minutes and 12 minutes, respectively.

###### 5.2.3. Optimization for All Ambulance Stations

Finally, we solve the concerned problem allowing all 17 potential bases to be optimized. In this context, the problem is feasible with the probation and the service reliability . After optimizing, can reach 0.71. The optimal objective value is shown in Table 6. Figure 7 shows the detailed optimal locations of the ambulance stations.

###### 5.2.4. Comparisons with Classical Models

We compare the model proposed in this paper with the classic DSM. This is because the model proposed in this paper is derived from classical DSM. The proposed model formulation considers two-stage rescue. We want to find out the difference between the classic DSM and our proposed model. However, in our model, we consider service reliability, which aims to find the optimal deployments of EMS vehicles to achieve maximized coverage of demand points. Each covered demand point is serviced by at least one ambulance in of the time. In order to examine the effect considering two-stage rescue, our model is slightly modified by removing the reliability constraint. Results are shown in Figures 8 and 9. The proposed model can achieve a better objective value than the DSM, with the primary coverage minutes. An interesting phenomenon is that when the number of stations exceeds eleven, our model can get full coverage. Besides, the coverage value obtained from the classic DSM cannot increase even if the number of stations exceeds nine. As for the range of primary coverage standard , our model can achieve full coverage in 8 minutes. However, the classic model gets full coverage with 12 minutes.

#### 6. Conclusion

This paper considers the optimization of the ambulance station location and ambulance allocation with the existing hospital locations to meet the patients’ injury types. After studying and analyzing the optimized ambulances locations for the Emergency Medical Services (EMS) in Songjiang, Shanghai, China, we concluded the ambulance base locations and the destination hospital locations affect the EMS performance significantly. Our proposed model can achieve remarkable results of coverage compared to the classic DSM in terms of used ambulance stations and the primary coverage standard . However, our results are based upon the real demands in only one-year data and our model is a deterministic model. It may be helpful in the strategic level but it lacks certain flexibility in the operational level. Also, we do not consider the cost of the ambulance station set up and the ambulance vehicle utilization.

Future research direction may include (i) the development of efficient algorithms for the large scale demand points and the potential ambulance station locations; (ii) the consideration of the EMS problem in a stochastic setting, such as considering stochastic emergency calls arrival and ambulance uncertainty.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Science Foundation of China (NSFC) under Grants 71571134, 71571135, and 71428002. The authors thank Teresa Tinta from University of Maryland, College Park, USA, for her help with the improvement of English writing.