key: cord-0733751-sxkuyf0e authors: Savaşer, Sinem Kınay; Kara, Bahar Yetis title: Mobile healthcare services in rural areas: an application with periodic location routing problem date: 2022-03-15 journal: OR Spectr DOI: 10.1007/s00291-022-00670-3 sha: 5da9b6c629a0c3c958d2cf30f05f3dc9d89294e5 doc_id: 733751 cord_uid: sxkuyf0e In this study, we focus on the delivery of mobile healthcare services in rural areas, where doctors visit remote villages which do not have a healthcare facility nearby. The aim is to increase the accessibility of healthcare services for such population centers. We aim to determine the village assignments of the doctors, their monthly visit schedules and base hospitals where they start and end their tours. We model this as a periodic location routing problem and use the policies of Ministry of Health of Turkey as a basis for our mathematical formulation. These policies include the essential components of mobile healthcare services, namely, continuity of care and determining evenly distributed periodic visits. We determine the visit schedules, i.e. routes, of doctors endogenously while satisfying these policies. We also develop a heuristic algorithm based on a cluster first-route second approach and solve larger instances more effectively. The computational experiments support that this solution methodology can effectively find optimal or near-optimal solutions and improve the computational times significantly. Availability of healthcare facilities, which may include multi-specialty hospitals or medical clinics, is one of the essential components of national health and welfare. The count and spatial distribution of these facilities have an impact on the widespread availability of healthcare services. However, convenient access to healthcare facilities in rural areas and by elderly patients remain to be a growing challenge. Traveling significant distances to reach the nearest healthcare facility to seek a simple medical treatment or to get routine checkups can be a real burden for these group of patients. The most effective solution to this globally recognized issue is considered to be providing mobile/home healthcare services, where healthcare providers visit patients instead of patients visiting healthcare facilities. There are many applications of mobile/home healthcare services throughout the world including the U.S., Canada, UK, France, Norway, Brazil, Egypt and India. It is observed that each country uses these services for different beneficiaries such as elderly, children, poor or patients living in rural areas (Sargutan 2006) . In the U.S., 12 million patients depend on some form of home care each year (NAHC 2020) . Similarly, 2 million Canadian patients, which is over 5% of their nationwide population, benefit from home care services annually (Statistics Canada 2020). Majority of the patients who are in need of these services are indicated as the elderly. The need for home healthcare services will continue to grow as the population projections show that we are moving towards an ageing world and 65+ population is expected to increase steadily within the next 50 years (US Census Bureau 2018) . In addition to these, recent developments show that home healthcare might have an increased importance during pandemics. The impairment of immune function of elderly people and likelihood of them to have co-morbid diseases (e.g. cardiovascular, respiratory, etc.) put them under the risk group. Home healthcare is expected to expose far less people to infectious diseases compared to their hospitalization. Recent COVID-19 spread has shown that it is essential to provide these patients the required support while they remain safely at home so that their public exposure at healthcare facilities is prevented. These services also help free up hospital beds and equipment during the pandemic for the patients who are in critical condition or in need of special medical equipment, e.g. respiratory ventilators. Besides the traditional home healthcare services, mobile practitioners/nurses visited vulnerable and high-risk group patients at home for COVID-19 testing and vaccination. At the beginning of the pandemic, testing via mobile units have been made available for those who are vulnerable and unable to get to an assessment center in Canada (Ontario Health 2020). In the U.S., Centers for Disease Control and Prevention published guidelines for mobile practitioners to vaccinate people who are homebound, i.e. who cannot leave home due to an illness (CDC 2021) . In Turkey, elderly who are above 85 have been granted the option to be vaccinated at their homes (CNNTurk 2021) . With similar practices all around the world, it is evident that mobile healthcare services could be essential for reasons other than traditional home healthcare during pandemics. Consequently, home healthcare operations are essential to support elderly, increase the access to healthcare in rural areas, and mitigate hospital bed and equipment shortages under pandemics. Therefore, it is crucial to assist healthcare decision makers in providing efficient home healthcare services. Motivated by this, in this work we consider the applications where family practitioners periodically visit remote rural communities (villages) which do not have a healthcare facility nearby. We develop a mathematical model that determines the monthly visit schedules of practitioners as well as their base hospitals while minimizing the travel time and satisfying problem specific requirements. Our study focuses on an application in Turkey where such a home healthcare system has been in use since 2010 in remote rural communities. There are two main policy requirements enforced by the Ministry of Health of Turkey in this application. First, the time between two consecutive visits to a specific village should always be the same, and we refer to this requirement as maintaining evenly distributed periodic visits. Second, the villages should be visited by the same practitioner, hence continuity of care should be satisfied (Official Journal 2010) . These type of periodic routing problems are already available in the operations research (OR) literature and they are mainly modelled with a predefined set of alternative schedules. The solution methodologies provided select the best visiting combination among many. However, unless all possible schedules are defined, which might be exponentially large, the solution of this approach might fail to find the global optimum. Another characteristic of these routing problems is that the servers are generally assigned to customers randomly, so that continuity of care is neglected. This is an important factor to consider in healthcare problems since it reduces fragmentation, thus, improves patient safety and quality of care. In connection to these points, the main contributions of this work can be listed as follows. We develop a new integer programming (IP) model which generates the schedules via its constraints and eliminates the need to define these schedules exogenously. These constraints are also tailored to meet the requirements of evenly distributed periodic visits. In addition to these, continuity of care is also considered in our modelling approach to increase the service quality of patients. This work fills a gap in the literature by providing a modelling approach that generates schedules, evenly distributes periodic visits and ensures continuity of care. We also develop a cluster first-route second based heuristic algorithm in order to obtain solutions in shorter computational times without compromising the solution qualities. We exploit the p-median formulation at the clustering phase of this algorithm and solve the routing subproblems to optimality by our proposed IP model. The remainder of this paper is organized as follows: in Sect. 2, we review the most relevant literature. In Sect. 3 we define the problem framework in detail and present our mathematical model. We explain the details of our heuristic algorithm in Sect. 4. Section 5 is dedicated to the discussion of the results of our solution methodologies. The study concludes with an overview of the work done and a discussion of adaptability to other possible application areas in this context. Home healthcare (HHC) is a well studied topic in the OR literature and it is introduced by Begur et al. (1997) and Cheng and Rich (1998) . The former study is motivated by a case study in the US and it develops a spatial decision support system to determine the schedules and routes of nurses. The latter work models the same problem as a multi-depot vehicle routing problem (VRP) with time windows and propose a two-phase heuristic algorithm that builds and improves tours at respective phases. Many other studies follow these two pioneering works and the majority of them focus on short-term or single period planning problems. As our study includes periodicity in the visits, we discuss the relevant works that consider multiple periods. For a detailed discussion on HHC routing and scheduling problems, we refer the readers to the review papers (Fikar and Hirsch 2017; Cissé et al. 2017) . In HHC problems, the most common issues addressed are determining the schedules and routes of the nurses. Steeg and Schröder (2008) is the first study that solves a multi-period HHC problem. The aim is to maximize the nurse-patient loyalty and the authors develop a hybrid approach that integrates adaptive large neighbourhood search (ALNS) and constraint programming (CP). Nickel et al. (2012) generates weekly schedules of the nurses with a two-stage algorithm that also benefits from CP. Grenouilleau et al. (2019) develops a heuristic based on set partitioning and ALNS frameworks. Braekers et al. (2016) extends the multi-period HHC problem by a bi-objective setting, where the travel cost and patient inconvenience are minimized. Cinar et al. (2019) introduces a priority parameter which increases exponentially by day for the unvisited patients. The authors simultaneously maximize the total priority of visited patients and minimize the travel time and solve the problem by an ALNS algorithm. Robust optimization is another OR tool that is used to address HHC problems. Cappanera et al. (2018) and Lin et al. (2018) introduce stochasticity to HHC problems by considering uncertainty in demand, e.g. last minute cancellations. Nasir and Dang (2018) considers patient and nursing staff selection according to the dynamic arrival and departure of patients. In addition to the classical HHC decisions, the authors also decide whether to accept a patient or hire an additional nurse. Demirbilek et al. (2019) also works on this dynamic acceptance and scheduling problem and the objective is to maximize the average number of daily visits of the nurses. A novel consideration of this work is to visit the patients at the same time of each week. There are a few studies which consider location decisions in addition to the classical decisions of HHC problems. Fard et al. (2018) is the first study to include location decision of the pharmacies that the nurses will be assigned to besides their routes. Fathollahi-Fard et al. (2019) extends this by considering emissions as a secondary objective and solving the problem for multiple periods with simulated annealing. As the routing component is the most prominent element of HHC problems, most of the developed models are adaptations of VRP. In almost all of the aforementioned studies, patients have time windows while they may report nurse preferences, and nurses are allocated to patients according to their qualifications. Except the last two studies, location decisions are not a part of the problem and continuity of care is overlooked. In the problem setting of our work, we focus on a more high-level problem where the routes of practitioners are determined on the village-level rather than patientlevel. As the practitioners are assumed to visit a group of residents at the same time slot at each village visit, this high-level problem setting does not include some of the main characteristics of the previous studies, e.g. patient time windows, nurse preferences, and nurse qualifications. In addition to this, our study differs by its periodicity of the visits (visit frequencies and rules) and continuity of care requirements. Thus, we believe Periodic VRP (PVRP) and Periodic Location Routing Problem (PLRP) are the most relevant modelling approaches to our problem setting, so we review this literature next. In classical PVRP, customers are visited multiple times during the planning period and the visit schedules are selected from a predefined set of alternatives. The aim of the PVRP is to find the allocation of customers to these predefined schedules such that each node is visited required number of times while minimizing the total cost. This problem is introduced to the literature by Beltrami and Bodin (1974) . The authors are motivated by a periodical municipal waste collection problem in New York City and prove that the problem is more difficult and complex than classical VRP. Russell and Igo (1979) proposes a heuristic approach which first selects a schedule for each node and then routes them individually. Christofides and Beasley (1984) provides the first mathematical formulation name the problem as The Period Routing Problem. This problem could not be solved optimally due to its complexity at that time. Once the problem is defined and proved to be NP-hard, heuristic algorithms dominated the literature in the succeeding years. Some of the articles address the problem by utilizing LP relaxation in their formulation (Tan and Beasley 1984; Chao et al. 1995) . IP based heuristics are also proposed, in which a part of the problem is solved via an IP and routes are generally constructed using Clarke and Wright algorithm (Russell and Gribbin 1991; Shih and Chang 2001; Gulczynski et al. 2011) . Tabu search (TS) is another commonly used algorithm to find higher quality solutions (Cordeau et al. 1997; Cordeau and Maischberger 2012; Liu et al. 2014) . There are also studies which propose metaheuristics that are based on TS (Archetti et al. 2018; Nair et al. 2018) . A column generation framework is used by multiple studies to solve the routing problem individually after assigning customers to the schedules by giving priority to the ones with the highest demand (Mourgaya and Vanderbeck 2007; Cacchiani et al. 2014) . Lagrangian relaxation, variable neighbourhood search (VNS) and genetic algorithm are the other approaches that are utilized for PVRP (Francis and Smilowitz 2006; Hemmelmayr et al. 2009; Carotenuto et al. 2015) . For a more detailed review of PVRP literature and used solution methodologies, we refer the readers to the surveys by Francis et al. (2007) and Campbell and Wilson (2014) . One of the common aspect of these studies is the absence of continuity of care and visiting rules. Rodriguez-Martin et al. (2019) is one of the few studies that implicitly incorporate continuity of care in PVRP so that the customers are served by the same driver at all visits. Another distinguishing factor of the previous studies is that they all define a set of alternative visiting combinations and their solution methodology chooses the best alternative among the ones in this set. However, it is time consuming and inefficient to define all such combinations, especially when the number of them increases exponentially. Hence, it is easy to see that unless all of them are defined exogenously, the solutions may be suboptimal. There are a few studies in the literature highlighting this issue under various applications of PVRP. Maya et al. (2012) works on a teaching assistant scheduling problem for disabled students in Netherlands. The authors take some visiting rules into account; however, do not impose continuity of care. An et al. (2012) determines the schedules of the nurses that are providing HHC in Korea. It considers continuity of care, but not include any visit rules, e.g. evenly distributed periodic visits. Archetti et al. (2017) defines the problem for a generic distribution plan where the quantity of delivered products at each time period is a decision variable. This study neither incorporates continuity of care, nor any visit rules. Thus, it can be claimed that none of the studies in the PVRP literature fully overlap with the requirements of our problem setting. The aim of the PLRP is exactly the same as PVRP, except it also finds the locations of the depots where the vehicles start and end their tours. The PLRP literature is not as broad as PVRP as it has not been studied until much recently. Since PLRP is even more complicated than PVRP, it is also NP-hard and the literature consists of heuristic approaches rather than exact solution algorithms. Note that, the PLRP setting is the one that is the closest to our problem framework. PLRP is introduced by Prodhon (2007) and a metaheuristic approach is proposed to solve it, in which locations and routes are determined, and improved by a local search. Prodhon and Prins (2008) proposes a memetic algorithm with population management for PLRP. It is observed that this methodology outperforms the previous iterative heuristic approach both in terms of solution quality and computational times. Prodhon (2009) and Prodhon (2011) develop an evolutionary local search (ELS) algorithm. Both of these studies seek to find the optimal periodic decisions by improving the assignments of the customers to the visiting combinations. Prodhon (2011) is the first study to develop a mathematical model with predefined schedules. It is indicated that the IP is capable of solving small instances and heuristics are essential to handle the large ones. Pirkweiser and Raidl (2010) addresses the same problem and develops a heuristic based on VNS. In the solution algorithm, the authors change depot locations and visit combinations iteratively and modify the daily routes by removing and reinserting the customers. Hemmelmayr (2015) utilizes a modified LNS algorithm which simply destroys a solution by removing customers and repairs it by adding them to another route or visit combination. It is observed that this approach improves the solution quality PLRP instances significantly. Koc (2016) introduces heterogeneous PLRP and its variant with time windows, and solves the large scale instances by a Unified-Adaptive LNS (U-ALNS) metaheuristic. The most recent study on this topic is by Hemmelmayr et al. (2017) . In this study, a similar approach to Hemmelmayr (2015) is followed and in each iteration of ALNS, the current solution is partly destroyed by a destroy operator and reconstructed by a repair operator. The classification of all PLRP studies can be found in Table 1 . This table points out that none of these studies take into account the three main aspects we consider in this study. Even though the characteristics of this problem are addressed separately in several studies in the PVRP literature, none of them considers creating the schedules with a mathematical formulation, satisfying visit rules as well as having continuity of care. Hence, it could be said that the dynamics of this problem generates a necessity for a novel solution approach. Before providing our IP model and solution methodology, we first discuss the problem framework and application specific requirements in the next section. We consider a problem framework where a healthcare decision maker seeks to find cost effective monthly visit schedules of practitioners (doctors) who are delivering mobile healthcare services in rural areas. The decision maker would also like to determine the base healthcare facilities (hospitals) where the doctors start and end their tours for reporting. There are additional application specific requirements such as evenly distributed periodic visits and continuity of care which are determined by the Ministry of Health of Turkey Official Journal (2010). According to the ministry guidelines, doctors are required to travel to the remote villages which do not have any healthcare facility nearby. Required monthly visit duration of a village depends on its population and for each duration, there are alternative visiting rules. Under our modelling framework, a working day consists of two half-day periods: 8am-12pm (morning) shift and 1pm-5pm (afternoon) shift. We then represent the monthly visit duration of each village as the number of half-days it has to be visited. This parameter will be referred as visit frequency hereafter. The relationship of population size and visit frequencies, as well as their alternative visit rules are shown in Table 2 . A small village with less than 100 residents has to be visited one half-day per month. For larger but still mid-sized villages, the visit frequency can increase to two, four or eight half-days in a month. The monthly frequency becomes 12 halfdays for the largest villages with more than 1,000 residents. Note that patients may also receive healthcare service without an appointment. Thus, doctors are required to be present at these villages even if there is no predetermined demand at a certain period. Selecting at least one visit rule among the alternatives ensures that visits are evenly distributed through a month to have a balanced schedule for the patients. The largest villages can only be visited in three consecutive half-days each week. A village with visit frequency of four (eight) is visited in one half-day each week (every 2.5 days) or in two consecutive half-days, which corresponds to a full working day, every 2 weeks (1 week). With a very similar logic, a village with visit frequency of two is either visited in one half-day every 2 weeks or in a day. Naturally, the smallest villages can be visited at any available time in a month. In addition to these rules, the villages have to be visited at the same slot in case they are visited in multiple weeks. For instance, a village with frequency of four visited on only Tuesday afternoon in the first week cannot be visited at any other time in the following 3 weeks. Similarly, a village with frequency of two can be visited on Friday morning of the fourth week if and only if it was visited on the same time of the second week. Besides reinforcing even distribution of the visits, this additional rule also ensures recurrent visits occur on a pattern so that the patients at a particular rural community can easily follow the arrival plan of the doctor. To demonstrate the frequency and visiting rules, we provide an example schedule with 8 villages, denoted by letters A-H. Assuming a work week consists of 5 days and a month is 4 weeks, there are 40 half-days in a month. Table 3 represents a schedule of a doctor that complies with the visiting rules: In this example, village A with frequency of 12 is visited in the first three consecutive half-days of each week. Villages B and C with frequency of 8 are scheduled according to the following alternative visiting rules: B is visited in every 2.5 days for a half-day throughout the month, and C is visited on Wednesdays every week. Among the villages D and E, which are the ones with frequency of 4, the former is scheduled to be visited on Thursday of weeks 1 and 3 whereas the latter is visited Week only on the Friday afternoons. Finally F is scheduled on Thursday of week 2, and G and H are visited in one half-day. Another essential requirement to satisfy is continuity of care. This ensures that each doctor visits the same villages during a month. Therefore, patient safety, familiarity and quality of care can be improved as the doctor knows the medical history of the patient. Consequently, symptoms and personalized treatments can be tracked efficiently in the following weeks. Finally, the base hospitals, where each doctor is going to start and end their weekly tours, have to be selected among the existing hospitals. We can summarize the unique characteristics of this problem as follows: • Visit frequencies depend on the population sizes. • There are alternative visit rules for each frequency level. • Services must be provided at the same slot each week. • Continuity of care is ensured. • Base hospitals have to be selected for each doctor. There are two assumptions we make while considering these requirements. First, we assume that the doctors stay at the villages they visit during the weekdays. This assumption is a reflection of the real-life application as the doctors receive a per diem by the Ministry of Health of Turkey, which is considered to be a significant incentive to work away from home throughout the week. Second, the start time of the afternoon shift is flexible. Since the doctors are staying at the villages after the afternoon shifts, these shifts do not have a strict end time. Therefore, if the travel time from a village to another is slightly longer, we assume that the afternoon shift can start and end later than desired. This is not expected to create significant delays since we will focus on a small region where the villages are not that far from each other. Considering these requirements and assumptions, we develop a mathematical model and a heuristic methodology that aim to determine the monthly schedules of the doctors and their base hospitals while minimizing the total travel time, which will be both cost and manpower efficient. In this section we introduce an IP formulation which determines the schedules of the doctors via its constraints. Before presenting the optimization model, we introduce the notation to be used hereafter: I Set of villages. It should be noted here that the capacity of each doctor is set to 40 time periods (half-days) assuming a 4 week horizon with five workdays. Time period 41 is defined to track the travel between the village visited on the last Friday afternoon and the base hospital. The decisions to be made can be represented by the following sets of binary variables: The following IP for the PLRP can now be proposed: The objective function (1) minimizes the total distance traveled by all doctors. In the first part, the distances are summed up for each trip from node n to m. In the middle part, distances between the villages where a doctor visits at the end of a week and at the beginning of the succeeding week are subtracted as doctors are routed between these points at times in T M in our modelling approach. In the final part, the distance between the node visited on Friday afternoon and the base hospital, and between the base hospital and the village that will be visited the following Monday morning are added up. By this way, total distance of doctors is calculated by considering their trips within a week as well as their weekly returns from/to the hospital. Constraint (2) guarantees that each doctor starts her tour at the beginning of the month. Each village is visited according to its visit frequency by constraint (3). Constraints (4) and (5) are the flow balance constraints. Unless a doctor is assigned to a village, she cannot be in that village at time period t due to constraint (6). Constraint assigns each village to exactly one doctor and ensures the continuity of care together with the constraint (6) which prevents visiting the villages with a different doctors. A doctor can be in or travel to at most one village at a certain time period by constraints (8) and (9). The working hour limitation of the doctors is satisfied via constraint (10) and it is guaranteed that they return to the hospitals as soon as their schedule is completed by constraints (11) and (12). Constraint (13) linearizes the k variable which can be defined as the multiplication of y and u. The location decisions are made via constraints (14)- (18). In this set of constraints, p hospitals are selected as bases, and every doctor is assigned to one. Any trip from a hospital is prevented if a doctor is not assigned to that hospital. The remainder of the formulation are the specialized constraint sets for each visit frequency and rule. For a village with visit frequency of 2, set of constraints (19)-(23) guarantee that once it is visited at time period t, then it should be visited either at the preceding or succeeding period (in case of two consecutive period visit), or 20 periods later (in case it is visited at the same half-day every 2 weeks). Note that, end of the weeks should be handled carefully to prevent the undesired schedules. For instance t = 10 and t = 11 cannot be accepted as two consecutive periods since they belong to different weeks. In order to prevent this type of issues, constraints (21) and (23) are defined separately. Similarly sets of constraints (24)-(28), (29)- (33) and (34)-(40) are defined to satisfy the visit rules of frequency levels of 4, 8 and 12, respectively. Note that for the villages with visit frequency of 4, bi-weekly schedules are constructed as the next 2 weeks of the month will be the repetition of the first two. Similarly, for the villages with frequency of 8 and 12 only weekly schedules are determined and they are replicated for the remainder of the month. Finally, the domains of the decision variables are given in constraint (41). Proposition 1 Even though there are two alternative visit rules for the villages with frequency of eight, they are visited in two consecutive half-days in an optimal solution as long as triangle inequality holds. Proof Assume to the contrary that there is an optimal solution in which a village with visit frequency of 8 is visited a half-day in every 5 periods. Consider an alternative solution which is obtained by moving the second visit of that village right after the first one in the first week, i.e. visiting this village in two consecutive periods while shifting the remaining villages by one period. It can be shown that this modified solution always provides a lower travel distance when triangle inequality holds. Since this village must be visited at the same time periods in the following 3 weeks as well, we can repeat this process and improve the objective function value further. Therefore, the initial solution cannot be optimal which creates a contradiction ◻. Next, we develop some valid inequalities to strengthen our IP formulation and decrease the computational times. These inequalities are logical derivations generated from the requirements of the problem and they can be represented as follows: Constraints (42) and (43) are implied by the capacity of the doctors. A doctor cannot be assigned to more than 3 villages with visit frequency of 12, and 5 villages with visit frequency of 8 and 12. The capacity restriction is additionally defined via u variables in constraint (44). Finally, the routes from a village are prevented for a doctor unless he is assigned to that village in constraint (45). Since different combinations of valid inequalities may have different effects on the solution times, we aim to find the most suitable one among 16 possible combinations. In order to do so, we perform extensive computational analysis on 15 test instances, which are detailed in Sect. 5.1. All possible combinations of these four valid inequalities are tested on a total of 240 configurations. For each configuration, we record the optimal solutions and the run times. We find the most effective combination(s) by evaluating their computational performance in terms of their run times. In order to determine the valid inequalities that generally perform well, we compare each combinations with the original model (without any valid inequalities). Additionally, we make instance-based comparisons by evaluating how each combination performs against the best one identified for that instance. The results of our analysis are summarized in Table 4 . Further details can also be found in Savaşer (2017). As a result of these extensive analyses, we decided to keep three combinations; only (44), (43)-(44) and all four inequalities. These three combinations have either the most number of instances with the top solution time or perform close to the fastest instance on average. In addition, they are also frequently better than the original model without valid inequalities. Since these three combinations provide desirable outcomes for the determined performance indicators, we continue to work with them. In the next step, these three combinations and the original model are tested on the 10 instances of the medium data set whose details are given in Sect. 5.1. It is observed that it gets harder to solve the mathematical model optimally when the data set size increases. Hence, we limit each experiment with a 2-h time bound. At the end of this duration, we record the value of the objective function and the remaining optimality gap for each setting. According to the computational analysis, we observe that the original model finds the best objective value in only one of the 10 instances. Combinations only (44), (43)- (44) and (42)-(43)-(44)-(45) find the minimum distance value in 2, 3 and 4 instances, respectively. We also observe that the last combination has the lowest average optimality gap and outperforms the original model in 6 instances out of 10. We refer the reader to Savaşer (2017) for further details on this data set analysis. As a result, we select this combination as the most suitable one for our PLRP model and the final version of the model denoted by P can be represented as follows: In this section, we develop a new heuristic based on a cluster first-route second approach in which we can exploit the hierarchical relationship between the assignment and routing decisions. In the first stage of the heuristic, we assign villages to doctors without violating working hour limitations. In the second stage, we determine the optimal schedule of each doctor individually based on their villages assignments which minimizes total traveled distance. There are some studies in the literature that use two-phased heuristic algorithms. Archetti et al. (2018) generates a feasible solution in the initial stage with a MILP and improve it by TS in the second phase. An et al. (2012) addresses the problem of nurse scheduling by considering only high-frequency patients at the first phase with a MIP formulation, and integrates the less-frequent ones into the solution in the following stage. Even though the idea of benefiting from an optimization program is similar, in this study we consider a cluster first-route second approach and have problem specific requirements that must be satisfied. In the first phase of our heuristic, we determine the clusters for each doctor by an IP formulation that takes the p-median model as a basis. This formulation determines the villages and the base hospital assigned to each doctor, which are within the working hour limitations, while minimizing the overall distance within the clusters. The notation we use in the IP is consistent with the original model. The additional parameters and decision variables are given below: C: number of clusters, i.e. number of doctors x j = 1, if village j ∈ I is selected as a cluster origin, 0, otherwise. y ij = 1, if village i ∈ I is assigned to cluster origin j ∈ I, 0, otherwise. The following IP for the clustering stage can now be proposed: In the objective function (46), we minimize the total weighted distance within all clusters. We prioritize larger villages by multiplying the distance with visit frequencies. We note here that alternative weighting schemes can be used in the objective function. For instance, one may choose to assign a weight of four to the villages with frequency of 12 and eight as these villages are consecutively visited in each week, and a weight of f i to the remaining villages. In constraint (47), the overall village set is divided into C clusters so that each doctor is assigned to a single one. Constraint (48) ensures that each village is assigned to an origin, i.e. every village receives healthcare services by a doctor. It should be noted here that the origins are dummy location points defined to use the p-median formulation and have no effect on the results of the heuristic approach. Constraint (49) satisfies the capacity limitations of a cluster (doctor). By constraints (50) and (51), a village cannot be assigned to a non-origin node and it is assigned to itself if it is an origin. p base hospitals are selected in constraint (52) and a base hospital is assigned to each cluster in constraint (53). Infeasible base hospital assignments are prevented by constraint (54). Finally, constraint (55) presents the domain restrictions of the decision variables. After determining doctors and the villages assigned to them, we use a simplified version of P in the routing phase of the heuristic. We adapt the formulation by first eliminating the doctor index d and location decisions as they are determined in the first phase. We also update the set I, distances and visit frequencies so that they only include the information of the villages which belong to that cluster. Therefore, we can claim that the computational challenges of our PLRP model are significantly simplified. With this adapted model, we determine the optimal schedule of each doctor separately and the objective values are summed up to obtain the overall distance traveled by all doctors. After finding an initial feasible solution to the problem, we investigate the ways of improving the overall distance value. We focus on improving the clustering stage of the heuristic with an iterative approach since the routing phase is solved to optimality by a mathematical model. In this approach, we aim to find the next best cluster at each iteration. In other words, we determine the clusters of villages that provide the minimum objective value which is greater than the previous iteration at each step. Then, we generate the schedules of the doctors and calculate the total distance value. At the end of the procedure, the smallest overall distance value is selected as the solution of the heuristic. In order to determine the next best cluster, we introduce an additional constraint to the clustering model, which can be represented as follows: represents the objective function value of the previous iteration and > 0 is a small enough number that prevents bypassing a feasible solution. With this constraint, we set a lower bound for the objective function of the clustering IP and force the model to find a new set of clusters. Note that, value is set to 0 in the first iteration to guarantee the feasibility of constraint. Despite this constraint, we can still obtain the same clusters with a different origins in the next iteration. Since the cluster origin is only a dummy decision in this process, these iterations cannot improve the solution. In order to prevent these repetitions, we add another constraint to the model: The set S represents the villages that are assigned to the first cluster, which includes village 1. This constraint tells that the villages in the first cluster cannot be all in the same cluster at the next iteration. We add this constraint to the formulation at each iteration and keep the ones that are generated in the previous steps. It should be noted here that, this constraint might skip the optimal solution as it is only stated for the first cluster. When we schedule more than 2 doctors, it is possible to find the optimal solution where the first cluster is the same but remaining ones are different. Using constraint (57) might bypass the optimal solution, whereas not using it might result in unnecessary repetitions. Therefore, we execute the heuristic with two different settings; the one with constraint (57) will be referred as Heuristic-1 and other without it will be referred as Heuristic-2. The representation of Heuristic-1 is provided in Algorithm 1. Notice that the representation of Heuristic-2 can be obtained by removing step 4. It should be noted that constraint (56) eliminates the possibility of finding an alternative cluster with the same objective function value. Although it is theoretically possible, it is highly unlikely to find such alternative clusters if the distance values are accurately estimated with many decimal places. If that is not the case, equation (56) may be replaced with (57) instead of using both of them simultaneously in Heuristic-1. In this section, we analyze the performance of our model P and heuristics via extensive computational studies. First, we explain the characteristics of the data sets we use. Next, we solve P with CPLEX and perform sensitivity analysis on the parameters of the model. Finally, we evaluate the performance of the heuristic approaches by comparing their solutions with P and each other. We provide the solution times of all three approaches and the differences between the optimization model and heuristic results. In this study, we use a real life data set from Burdur, which is among the first cities that the Ministry of Health provided mobile healthcare services in Turkey. The data covers three main municipalities of the city which have 50 residential points in total. Three different sized sets are formed as small, medium and large data sets. These are generated by taking the first 15, 30 and all 50 nodes of the main data set; where 12, 26, 45 of the nodes are the villages and the remaining ones are alternative base hospital locations, respectively. The details of the data sets can be found in Kara and Savaser (2018) . We gather the exact coordinates of these villages from Google Maps and calculate the Euclidean distances between each pair. In addition, we determine the visit frequencies according to the population size of each point. However, we also generate alternative artificial frequencies to perform detailed computational analysis. We set the capacity of each doctor to 40 periods. Finally, we vary the number of base hospitals according to the number of doctors required for each instance. These result in 30, 10 and 10 instances for small, medium and large sized samples, respectively. In this part, we first focus on solving P using an off-the-shelf solver, CPLEX, and analyze the solutions. We perform sensitivity analysis in order to understand the changes in the solution times of P when a commercial solver is used. To do so, we vary the number of doctors (size of set D) and base hospitals to be selected (p). We perform additional sensitivity analysis by changing the distribution of visit frequencies (f). We solve P using CPLEX with all three data sets. For the small data set, we find the optimal solutions in reasonable times. Even though we try to solve the medium and large data sets with 2-h and 4-h time limits, respectively, we are unable to reach to the optimal solution due to the size of the problem in either case. Setting longer time-limits does not help as the solver terminates with an out-of-memory error. We present the findings on the small data set in Table 5 . According to this table, the computational times show significant variations from each other in 30 instances. The mean of the solution times is 2,811 s which is slightly above 45 min. The standard deviation of the instances is 2,984 s. The 90% confidence interval of the solution time of this data set is within 1,915 and 3,707 s, i.e. 32 and 62 min. Table 5 is going to be used in Sect. 5.3 while comparing the heuristic results with the optimal solutions. The detailed instance specifications are provided in Table 8 , and the correlation between different solution times and instance characteristics are discussed thoroughly thereafter. We then analyze the effects of the parameters on the solution times. We triplicate each instance by changing the order of the visit frequencies to increase the size of our data set. Therefore, the number of villages with each frequency level remains the same; however, the distribution of frequencies among all villages change in the other two variants of that instance. The number of doctors required to visit all villages varies according to the frequency levels. More villages with frequency of 12 directly increases the required number of doctors. On the other hand, scheduling less doctors can be enough when the majority of the villages are smaller with frequency of 1 or 2. In our data set, a single doctor is scheduled in 6 of the instances. There are 18 and 6 other instances that require to schedule 2 and 3 doctors, respectively. After solving P for all these instances, we determine the minimum, average and the maximum solution times according to the varying number of doctors. These results can be found in Table 6 . It can be observed from the results that increasing the number of doctors makes the problem more complex as expected, thus, the solution times increase dramatically. While the schedule of a single doctor can be determined to optimality in less than a minute on average, the schedules of multiple doctors can be generated in Next, we systematically vary the number of base hospitals (p) over the 30 instances of the small data set. In fact, certain instances are duplicates of each other with different p values. For instance, instances 3 and 4 have the same frequency distribution and three doctors; however, instance 3 selects two base hospitals, whereas instance 4 selects three. There are 18 instances that pick one, 10 instances that pick two and 2 instances which pick three base hospitals. In order to make more reliable comparisons, we categorize the instances according to the number of doctors and then perform the analyses on the number of base hospitals. The respective solution time analysis can be found in Table 7 . In Table 7 , it is observed that selecting less base hospitals results mostly in shorter solution times. The instances that have only one doctor cannot be compared with other instances as they can only select one base hospital. However, they provide the shortest solution times on average, which is parallel to the outcomes of the previous analysis. When we schedule two doctors, p = 1 provides a lower average solution time than p = 2 . On the other hand, there could be exceptions as the maximum solution time of the former setting is higher than the latter one. We obtain similar outcomes when we route three doctors. When there is only one base hospital for all doctors, the model can be solved in around 1.5 h on average. Approximately 18 additional minutes are required to solve the setting which selects one more base hospitals for three doctors. When p = 3 , the average solution time becomes almost 2.5 h. There are also exceptions in scheduling three doctors that are observed in minimum solution times; however, the computational times increase with the increasing number of base hospitals in majority of the cases. We then analyze the distribution of the workload between doctors. To do so, we first assume a balanced allocation is obtained when the monthly workload difference between any two doctors is less than or equal to 6 half-days. When there are minimum number of doctors required, i.e. �D� = � ∑ i∈I f i 40 � , and total workload of the 3,765.36, 13,523.98] system is at least 90% of overall doctor capacity, 90% of our solutions provide a balanced allocation. On the other hand, if there are more than enough doctors to schedule, the results may become unbalanced. For instance, when there are two doctors for a total workload of 40, a solution that allocates 30 half-days to one doctor and 10 half-days to another may be optimal. This is because we do not model workload fairness between doctors explicitly in our approach; however, this would be an interesting future research direction. According to this analysis, we recommend using our model when the number of doctors is sufficient but not excessive, and the monthly workloads would plausibly be fair. As it is mentioned above, certain instances are duplicates of each other with different p values. We also compare these pairs with each other in order to observe the effects of varying number of base hospitals for the exact same setting. These comparisons can be found in Fig. 1 . There are eight instance pairs that are duplicates of each other. Among these instance pairs, we find out that increasing the number of base hospitals (p) by one unit strictly increases the solution times, either slightly or dramatically. For example, instance 19 requires only 7% more time than instance 18, whereas the solution time of instance 14 is almost 15 times higher than the computational time of instance 13. Therefore, we can claim that higher p values for the same instance increases the complexity of the problem and results in higher solution times. We also investigate the allocation and routing decisions in these duplicate instances when p is increased. Among the eight instance pairs, village allocations to doctors do not change in five of those when there is an extra base hospital. In two of those instances, we observe only minor allocation changes, e.g. the allocation of two villages are swapped between doctors. On the other hand, we observe that base hospital and doctor allocations may change when an additional hospital is open. Particularly, this allocation change is more likely to happen if the additional base hospital is located closer to the villages a doctor visits. While solving Ins10 - large scale instances, this observation may be leveraged to determine good solutions without solving the entire problem again. Instead, a doctor may be allocated to the new base hospital, especially if the villages are closer to that location, and the route of that doctor can be optimized quickly. We then evaluate the frequency distributions by investigating the impacts of having a dominant frequency According to the 90 settings of the small data set, average solution times of 3 variants of 30 instances are provided in Table 8 . Before interpreting the solutions, we categorize the instances according to the visit frequency that is observed at least six times, so that this frequency level constructs the majority. According to this, the instances that have majority of frequency of 12 are 1-6. The common characteristic of these instances is that they have to schedule three doctors. Therefore, aligning with our previous observations, the solution times of these instances are the highest with the exception of instance 2. On the contrary, the instances that contain mostly frequency of two and one, which are instances 25-30, require only 1 doctor. Hence, they are solved in less than 100 s on average. The instances with majority of frequency of eight are 15-19 and four are instances 20-24. We observe that these instances mostly provide results within 2,000 and 4,500 s, which indicates that they cannot be solved easily. On the other hand, the remaining instances 7-14 do not have any dominating visit frequency, i.e. the frequencies are balanced among all villages. When these instances are analyzed, we find out that they can be solved in shorter times compared to others with the exception of instance 14. Moreover, among these 8 instances, {7, 8, 10, 13} select only one base hospital and they are solved under 5 min. Overall, we can say that having villages with visit frequency of one and two solves the model faster, whereas dominance of frequency of 12 increases the solution times. This situation can be easily associated with the varying number of doctors. Besides that, we can claim that having evenly distributed frequencies and selecting less base hospitals result in shorter computational times and eases the complexity of problem significantly. In order to derive some insights about features of an optimal solution, we finally analyze which of the alternative visiting rules are selected in optimal solutions for different visit frequencies. In accordance with Proposition 1, villages with frequency of eight are always visited in two consecutive half-days every week. In 30 instances we look at, there are 82 and 55 villages in total with frequency of four and two, respectively. 78% and 87% of the time, these villages are visited in two consecutive half-days in the same week to minimize the total travel time. We observe that these villages are generally visited in one half-day per week (or 2 weeks) when the remainder of the villages that are assigned to the same doctor have higher frequencies. For instance, assume a doctor is allocated to three villages with frequency of 12 and a single village with frequency of four. In this case, there are no feasible solutions where the village with frequency of four is visited consecutively in a week, hence non-consecutive visits are necessary. As a result, we observe that the villages with alternative visiting rules are most likely visited in consecutive half-days as long as such a visit schedule is feasible. This outcome may be utilized by the decision makers to determine good solutions for large scale instances when using optimization models are time-intensive, impractical or impossible. The main goal of the heuristics is to determine solutions in shorter times than the mathematical model without significant solution quality compromises. With this purpose, we first solve 24 instances of the small data set with both heuristic settings and compare their outcomes with the optimal results of P . We set the number of iterations to 20 as it provides high quality solutions in short amount of time. Note that solving instances 25-30 of this data set cannot provide any insights about the heuristic performances. This is because they schedule only one doctor and the first iteration of the heuristic will be exactly the same as solving this instance with the original model. The results of Heuristic-1 and Heuristic-2 are presented in Table 9 . For each variation, we provide their objective value, solution time, and optimality gap. We also give the iteration number that provides the best result for each instance under the 'Iteration' column. To compare the computational performance of the heuristics with the solver, we provide the difference between their solution times, where negative values show by how much the heuristics are faster. We observe that both heuristic variations provide significantly better computational times than the mathematical model at each instance, where they are approximately 56 min faster than the solver on average. Heuristic-1 provides a solution between 1 and 4 min with the average of 132 s and Heuristic-2 solves the problem within 1.5 and 5 min while its average is slightly less than 3 min. Another observation is related to the solution qualities. Among the 24 instances, Heuristic-1 is able to find the optimal solution in 18 of them, whereas Heuristic-2 can find it in 16 instances. The average optimality gaps are only 0.79% for the first variation and 1.12% for the second one. Thus, it can be claimed that both approaches are eligible to find high-quality solutions in the majority of the instances and perform quite well on the small data set. When we compare the two variations, Heuristic-1 seems to perform better in terms of solution time and quality. However, there are still some instances where Heuristic-2 outperforms the first one. In order to illustrate both cases, first consider instance 14, where Heuristic-1 finds the optimal solution; however, Heuristic-2 fails to do so and finds a solution that falls 6.88% behind the optimal one. The main reasoning behind this is the additional constraint (57). Since repetition of the clusters are not allowed in this approach, Heuristic-1 determines different cluster sets and finds the best at the 16th iteration. However, the same clusters are observed multiple times in Heuristic-2 and the optimal clusters cannot be found in 20 iterations. On the other hand, this heuristic finds the optimal solution in instance 19 but Heuristic-1 cannot. As it is discussed in Sect. 4, constraint (57) eliminates the optimal clustering structure and finds a solution which has an 1.32% optimality gap. Therefore, we can see that both approaches have their own drawbacks in terms of finding the optimal solution. In order to provide more managerial insights, we map and compare the routes of doctors based on the results of the mathematical model and heuristic approach. We consider instance 16 and depict the routes of two doctors over 4 weeks in Fig. 2 . With a 2.58% gap from the optimal solution, heuristic assigns very similar villages to each doctor. The only difference is that second doctor (indicated with blue) is responsible for village 5 instead of 11 and 12. Other than that, it could be observed that the structure of the routes are similar to each other, indicating that heuristic solutions are still practically applicable. We proceed by testing the heuristic on the medium data set. We observe that the best results are found within the first 10 iterations of the small data set in 88% of the instances. Since the schedules are determined with the mathematical model, we decrease the number of iterations to 10 to keep the solution times reasonable. We record the same values for this data set; however, instead of the optimal objective value of P , we provide the value obtained at the end of 2 h. These results can be found in Table 10 . The instances can be solved in 171 s on average with Heuristic-1. It takes around 40 more seconds for Heuristic-2 to solve the same instances. In addition to this, both heuristics improve the solutions that the model finds in 2 h in most of the instances. We observe 5.73% and 6.22% improvements on the solutions on average by the first and second variations, respectively. Both approaches provide the same result for the first four instances, and among the remaining six, each setting outperforms the other in half of them. Similar to the small data set results, there are certain cases where one methodology provides better results than the other, hence, we cannot determine a strict domination between the two. Yet, we can claim that the heuristic approaches perform well and they provide better solutions in 3-4 min than the solver generates in 2 h. Finally, we solve the large data set. We provide the results of the mathematical model after 4 h as well as the heuristic distances, solution times and improvements over the solver solutions. The outcomes can be found in Table 11 . As the size of the instances grow, the quality of solutions decrease for the mathematical model. Therefore, we observe considerable improvements on the distance values in both Heuristic-1 and Heuristic-2 solutions. The average improvement percentages of respective methodologies are 62.6% and 62.5%, which do not differentiate them. As it is observed in both small and medium data set analyses, Heuristic-1 provides the results in slightly shorter times, namely in 6 min 21 s on average, whereas Heuristic-2's average computational time corresponds to 6 min 43 s. Among the 10 instances of the large data set, Heuristic-1 provides better solutions in 6 of them while Heuristic-2 outperforms the former approach in 3 instances. In the end, we cannot claim that one of the heuristic variations is significantly better than the other. Depending on the instance, either methodology can generate a solution that the other one skips. Nevertheless, we can still see that Heuristic-1 finds better solutions in slightly more cases. In addition to that, the solution times of this approach turn out to be better in almost all of the instances. Even though we develop Heuristic-2 as we suspect that initial setting would skip solutions that might be optimal, we find out that clustering stage performs well enough to find qualified solutions in shorter times in Heuristic-1. In this study, we focus on the mobile healthcare problem in rural areas, where practitioners travel between the villages to provide services and the aim is to increase the accessibility of these services in the underdeveloped communities. We introduce a model that considers essential requirements such as continuity of care and evenly distributed visits. We use the policies of Ministry of Health of Turkey as a basis for our mathematical formulation. Yet, our modelling approach can easily be adapted to similar policies of other countries. We present the first study in home healthcare and PLRP literature that determines the schedules of the doctors endogenously, ensures continuity of care and evenly distributed visits simultaneously. As determining the routes in the model increases the size and complexity of the problem, we also develop a heuristic based on a cluster first-route second approach to solve larger instances more effectively. In the initial phase of this solution methodology, we allocate villages to the doctors by an IP that is based on p-median formulation and we route the doctors with the IP model that we develop in the second phase. In order to decrease the size of the solution space and improve the computational performance, we propose valid inequalities and test their computational efficacy. We observe that the small sized sets can be solved in reasonable times using a commercial solver. Computational experiments and sensitivity analysis reveal important managerial insights about the correlation of instance specifics and optimal solutions. We find out that scheduling more doctors and selecting more base hospitals increase the solution times, whereas having more balanced visit frequency distributions helps finding optimal solutions in shorter times. As the IP model do not perform well when the sizes of the instances increase, we test and evaluate the performance of the heuristic. We find out that they reduce the computational times drastically and generate optimal or near optimal solutions in the small data set. In the medium and large sets, we observe that the heuristic improves the objective values of the mathematical model that are found in certain time limits in all instances and we obtained these improved solutions in much shorter times. Therefore, we can claim that we develop a fast and efficient heuristic for our problem. From a managerial point of view, the importance of monitoring the health of a community has been increasing since the COVID-19 pandemic. It is crucial to provide systematic healthcare services, especially for those who do not have easy access. In addition, ensuring continuity of care with dedicated patient-doctor assignments and following up status of patients are vital to improve the overall health of the population. Mobile services may also play a significant role during the pandemics in providing regular testing and vaccination efforts. Therefore, we develop an efficient algorithm as a tool for the decision makers so that they can rapidly find cost-effective schedules for the assigned doctors who are frequently visiting the rural communities. Even though we focus on delivering mobile healthcare services in rural areas, there are many other applications of mobile services under healthcare as well as different areas. For instance, nurses provide home healthcare to the elderly quite often, and specialists (e.g. cardiologists, dermatologists) visit patients that require continuous care and chronic treatment. We have also seen that the home healthcare becomes more important under unprecedented times such as a pandemic. Education or sanitation services could also be provided with periodic mobile services. Recently, as a result of the arising need through the world, doctors, teachers and sanitation engineers started to travel through the refugee camps to provide respective services. These problem types can be addressed by extending our work based on the problem-specific requirements of these applications. Another extension of our study is to take workload fairness among doctors into account. Although considering fairness among patients is a priority in the healthcare logistics applications, an interesting future research direction could be improving the quality and practicality of the schedules by incorporating workload fairness among doctors. Scheduling healthcare services in a home healthcare system The flexible periodic vehicle routing problem A two-phase solution algorithm for the flexible periodic vehicle routing problem An integrated spatial dss for scheduling and routing home-healthcare nurses Networks and vehicle routing for municipal waste collection A bi-objective home care scheduling problem: analyzing the trade-off between costs and client inconvenience A set-covering based heuristic algorithm for the periodic vehicle routing problem Forty years of periodic vehicle routing Demand uncertainty in robust home care optimization Periodic capacitated vehicle routing for retail distribution of fuel oils Vaccinating homebound persons with covid-19 vaccine An improved heuristic for the period vehicle routing problem A home health care routing and scheduling problem. Technical report Christofides N, Beasley J (1984) The period routing problem Or problems related to home health care: a review of relevant routing and scheduling problems Citizens over 85 will be vaccinated at home today A parallel iterated tabu search heuristic for vehicle routing problems A tabu search heuristic for periodic and multi-depot vehicle routing problems Dynamically accepting and scheduling patients for home healthcare A location-allocation-routing model for a home health care supply chain problem. World Academy of Science, Engineering and Technology A green home health care supply chain: new modified simulated annealing algorithms Home health care routing and scheduling: a review Modeling techniques for periodic vehicle routing problems The period vehicle routing problem with service choice A set partitioning heuristic for the home health care routing and scheduling problem The period vehicle routing problem: new heuristics and realworld variants Sequential and parallel large neighborhood search algorithms for the periodic location routing problem A variable neighbourhood search heuristic for periodic routing problems A periodic location routing problem for collaborative recycling Periodic location routing problem: an application of mobile health services in rural areas A unified-adaptive large neighborhood search metaheuristic for periodic location routing problems Jointly rostering, routing, and rerostering for home health care services: a harmony search approach with genetic, saturation, inheritance, and immigrant schemes Hybridization of tabu search with feasible and infeasible local searches for periodic home health care logistics A metaheuristic for a teaching assistant assignment-routing problem Column generation based heuristic for tactical planning in multiperiod vehicle routing Scheduling and routing models for food rescue and delivery operations Solving a more flexible home health care scheduling and routing problem with joint patient and nursing staff selection Mid-term and short-term planning support for home health care services Regulations of Family Medicine Practice Progress made on covid-19 specimen collection and lab-analysis Multilevel variable neighborhood search for periodic routing problems A metaheuristic for the periodic location-routing problem An elsxpath relinking hybrid for the periodic location-routing problem A hybrid evolutionary algorithm for the periodic location-routing problem A memetic algorithm with population management (ma|pm) for the periodic location-routing problem The periodic vehicle routing problem with driver consistency A multiphase approach to the period routing problem An assignment routing problem In: Sargutan A (ed) Comparative healthcare systems (in Turkish) Periodic location routing problem: an application of mobile health services in rural areas Receiving care at home A hybrid approach to solve the periodic home health care problem A heuristic algorithm for the period vehicle routing problem Older people projected to outnumber children for first time in U Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations The authors would like to thank the associate editor and two anonymous reviewers for their comments, insightful suggestions and careful reading of the manuscript. The authors declare that they have no conflict of interest.