Biodiversity constraint indicator establishment and its optimization for urban growth: framework and application

Urbanization causes tremendous pressure on biodiversity and ecosystems at the global scale. China is among the countries undergoing the fastest urban expansion. For a long time, ecological environment protection has not been a priority in China’s urban planning process. Current urban growth optimization research has some limitations regarding the selection of more scientific ecological constraint indicators and the interaction between urban expansion and ecological factors. This paper at first aimed to establish a reasonable comprehensive biodiversity constraint indicator based on the indicators of net primary productivity, habitat connectivity and habitat quality, and then conducted a case study in Beijing and compared biodiversity loss and urban growth patterns under different developing situations. The integrated valuation of ecosystem services and trade-offs model and GIS-related methods were used to obtain biodiversity and ecological spatial distribution layers. Then an ecological priority-oriented urban growth simulation method was proposed to search for minimum biodiversity loss. The results showed that the important biodiversity security areas were mostly distributed in the western part of the study area and that the ecological degradation in 2000 had a radial pattern and was well in line with the urban construction and ring road distribution patterns. Meanwhile, biodiversity loss with the biodiversity constraint was much less than actual urban growth in 2000–2010. Under the guidance of ecological optimization, urban growth in the research results reflects decentralized and multi-center spatial development characteristics. This type of urban growth not only provides a new model for breaking the inertia of urban sprawl but also proposes ‘biodiversity security’ as an applicable regulatory tool for urban planning and space governance reforming.

. Location of the study area. Within the study area boundary, 1 is the district of Haidian; 2 is Shijingshan; 3 is Fengtai; 4 is Dongcheng; 5 is Xicheng; 6 is Chaoyang.

Introduction
City growth and the urbanization process are strong linked with biodiversity and ecological processes [1]. Although urban land cover is a relatively small fraction of the Earth's surface, urban areas drive global environmental change [2]. Currently, urban expansion has an increasingly intensive impact on biodiversity and ecosystem services (ES) [3][4][5][6].
Today, urban areas globally are expanding on average twice as fast as their populations [7,8]. Nearly half of the increase in high-probability urban expansion globally is forecasted to occur in Asia, with China and India accounting for 55% of the increasing [9]. In urban planning, biodiversity and ecological factors are less important than other considerations [10][11][12]. For a long time, ecological environmental protection has not been a priority in China's urban planning process [11]. Synthesizing biodiversity and ecological factors with city growth and urbanization is essential for sustainable development. Developing an analysis of the interaction between urban expansion and the ecological environment has been a hotspot in geography and ecosystem studies [13].
Currently, many studies analyze the relationship between urban expansion and the ecological environment. Most focus on two aspects: the urban expansion impact on the ecological environment and the ecosystem restriction on urban expansion. Regarding the former, some studies concentrate on the direct habitat problem, such as natural habitat loss and degradation [9,12,14,15] while others use an ecological index to indirectly represent habitat loss and degradation, such as the landscape fragment index [15,16] or the habitat pressure index [17]. Recently, ES are increasingly used to reflect urban expansion's influence on ecology [16,[18][19][20]. Those studies select ES types according to their importance and the data collection situation in their study areas. Other studies select some typical species to show the ecology degradation information [9,12,21].
Many studies have also examined the ecological factors constraining urban expansion [22][23][24][25][26]. These studies mostly build landscape security patterns as the constraint factors impacting urban expansion, from the perspective of the ecosystem structure and its functional stability [27]. The framework of the built landscape security pattern has 3 steps: firstly, it defines the core ecological area, such as the protection area, the nature park, and the water resource area [28,29]. Secondly, it uses some models (such as minimal cumulative resistance) to calculate the block degree of all the land use types [25,30]. Thirdly, using the predicted amount of construction land to fill the less-valuable ecological areas in a bottom-up fashion, the urban expansion space can be obtained [28]. Urban expansion models such as cellular automata (CA) sometimes consider ecological constraint factors (e.g. [18,20,31]). However, these models normally simulate urban expansion according to some socioeconomic factors; they take ecological indicators as the block area for urban expansion but do not consider parameters attending the urban expansion.
The current research is limited in its selection of more reasonable ecological indicators. Those ecological indicators, especially the ES indicator, differs distinctly in most studies (e.g. [18-20, 32, 33]), which makes cross-regional application and comparisons difficult. Moreover, representatives of these indicators are uncertain and unreasonable. For example, the food supply is not a core ES in the urban ecosystem because normally it can be obtained from commercial sources [26]. However, the food supply is often included in current studies [19,20,32] as are other nonimportant urban ES. Meanwhile, some important services, such as cultural services, are usually ignored [20]. Species indicators may not be adequate or sufficiently accurate to represent the region's biodiversity [34] because whether some special species can represent the biodiversity of one large region is uncertain and data on species distributions are usually scarce (Christian et al 2015). Another problem with the related studies is that they normally do not pay much attention to the interaction of urban expansion and ecological factors. Some studies analyze the linear correction between the urban expansion result and the ecological indicators result [16,20], but they do not reflect the interaction between the urban expansion process and the ecological process. Even the urban expansion model hardly considers urban growth impacts on ecological processes [18,31], nor does it join the ecological process with the synthesis calculation in the urban expansion model; therefore, most current models cannot detect urban expansion from the perspective of ecological mechanisms and processes [18,20,35].
The current framework of the built landscape security pattern as the biodiversity and ecological constraint factors for urban expansion has some limitations as well. Under that framework, studies usually tend to pay more attention to the natural attributes of habitat but ignore the relationship between the habitat and surrounding environment [36]. Furthermore, studies usual build the block area for urban expansion and identify its value based on some coarse and subjective judgments. Finally, urban expansion is a complex socioeconomic and ecological phenomenon [37]; the expansion process is obviously accompanied by ecological environmental dynamic change. However, currently the block area in that framework is simply a static constraint for expansion that does not consider the interaction between urban expansion and the ecosystem.
Considering the limitations of current studies, this paper has the following aims: (1) establish reasonable biodiversity constraint indicators for urban expansion and analyze the foundation biodiversity security pattern of the study area. (2) Compare biodiversity loss characteristics under different urban expansion patterns, with or without the biodiversity constraint.

Study area
Beijing is the capital of China, located in northeast Asia. Since the reform and opening-up, Beijing has experienced tremendous urban expansion, along with rapid economic development and substantial population growth [20,38]. Since then, the urban land area in Beijing has increased by 4.5 times, the population has approximately doubled and the gross domestic product has increased approximately 14 times [20].
Beijing has four types of functional areas: the city center, expansion areas, new development areas, and biological conservation areas [19,39]. We choose the city center and expanding areas as our study area (figure 1), as these areas are the main urban growth areas. The administrative districts of the study area are Dongcheng, Xicheng, Shijingshan, Fengtai, Haidian, Daxing, Tongzhou, Chaoyang and Shunyi.

Methods
Establishment of a biodiversity constraint indicator Biodiversity and ES are intrinsically linked [40]. Biodiversity, including habitat diversity, contributes to human wellbeing by ensuring ecosystem functionality and resilience [41]. Biodiversity is also an important determinant of ecosystem processes and contributes to the provision of many ES [42]. From the landscape-scale perspective, high-level biodiversity can be defined from 3 aspects [27,30]. First, biodiversity can support high-level ES. Second, biodiversity can reflect the connectivity of habitat patches. Third, biodiversity can prevent habitat degradation. Specifically, energy is the basis of survival for all kinds of living beings, and net primary productivity (NPP) is an appropriate index for representing the energy produced by vegetation. Habitats with high productivity have more available sources and sustain more individual populations, increasing the survival probability of species [43][44][45][46][47]. Habitat connectivity is crucial for maintaining ecological processes in good condition. Habitat connectivity is also a key factor in maintaining the stability of biodiversity and ecosystems [48,49]. Habitat degradation determines habitat quality, and habitat quality depends on land use/cover structure and the level of intensity in the surrounding area [50].
In this paper, we established an integrated biodiversity constraint indicator at the landscape scale by integrating NPP, habitat connectivity and habitat quality (Formula1), in which each aspect was considered equally important because they are all crucial for biodiversity protection and ecological stability. The final value of this index was standardized into a range of 0-1 by the method of min-max normalization, with higher values representing higher biodiversity levels.
= + + ( ) BI NPP CON HQ. 1 In this formula, the BI (biodiversity constraint indicator) represents the biodiversity level at the landscape scale. We calculated the monthly average NPP of Beijing in 2000 as the NPP biodiversity constraint indicator. Then, we used the method of min-max normalization standardized the NPP value over the range of 0-1. CON represents the connectivity of each habitat patch, it was original ranged from 0 to 1. HQ refers to the habitat quality, which was obtained from the integrated valuation of ES and trade-offs (InVEST) model and ranged from 0 to 1.
Any area or region with exceptionally high biodiversity at the ecosystem, species and genetic levels can be defined as a biodiversity hotspot [51]. Biodiversity hotspot analysis was used to obtain the biodiversity security pattern of the study area based on the distribution of NPP, habitat connectivity and habitat quality.

NPP calculation
The net primary production of vegetation is the production of dry organic material of vegetation in set time and area. NPP can reflect the status and the flexibility of the ecosystem and it is the basic condition of survival for all kinds of living beings in an ecosystem [52]. It can be calculated by Carnegie-Ames-Stanford Approach model [53,54]. The calculation depends on the variant of vegetation's absorption of photosynthetically active radiation (APAR) and light utilization efficiency (ε), It can be expressed in equation (2): In this equation, NPP(x, t) means net primary production of vegetation on location x during time t. APAR (x, t) means vegetation's absorption of photosynthetically active radiation on location x during time t. ε(x, t) means conversion efficiency to solar radiation on location x during time t. ε(x, t) can be expressed as follows: In this equation, e max means the max conversion efficiency to solar radiation; T(x, t) means the temperature factor on location x during time t. W(x, t) is the humidness factor on location x during time t. Parameter in NPP calculation was set by the project of 'China's national ecosystem remote sensing survey of the year 2000-2010' [55].

Habitat connectivity calculation
Habitat connectivity is the degree of the landscape facilitates the movement of species and other ecological flows [49]. Saura and Pascual-Hortal developed the connectivity index 'probability of connectivity' (PC) to represent habitat connectivity conditions at the landscape scale [56,57]. PC is defined as the probability that two points randomly placed within the landscape fall into habitat areas that are reachable from each other. The value of PC is in the range of 0-1, a high value means high probability. The computation of the PC index benefits from the application of graph-based algorithms for determining the maximum probability paths [56], and it can be shown by the following formula: where n is the patch amount at the landscape scale and a i and a j are the areas of patches i and j, respectively. A L indicates the total area of the study area. * P ij refers to the expanding probability between patch i and j. In this study, PC is used to assess the connectivity ability of each habitat patch. Conefor Inputs for ArcGIS 10.0 (software package) and Cone for Sensinode 2.6 (software package) were integrated to calculate the PC index. In our study, we considered the area which PC value above 0.5 as the good habitat connection area.

Habitat quality assessment by the InVEST model
Habitat quality refers to the ability of the environment to provide conditions appropriate for individual and population persistence [58]. Habitats in a favorable conservation status host a richer biodiversity than those in an unfavorable conservation status [59]. In this study, we used the InVEST model (3.6.0, [50]) to evaluate the HQ of the urban area in 2000. This approach is particularly suitable for studying HQ at a large scale and has been proven to be effective in characterizing biodiversity [60][61][62][63]. HQ strongly depends on the suitability of each land use and land cover (LULC) class to provide habitat for wildlife. The mechanism of HQ evaluation in InVEST considers the external stress (threats) impact on wildlife habitat. Threats factors include: intensity of each threat, distance between the habitat and the stressor, the associated stress pathways (both linear and nonlinear), and habitat sensitivity to each threat [50]. Threats and habitat are represented by LULC data. In our study, threats are represented by the urban area, the village area, the construction land outside the urban and village areas, the agriculture area, the railway and all kinds of roads. Habitat is represented by forest land, grassland, farmland, wetland, and some other unused land for wildlife. The HQ score can be expressed as in equation (4): where Q xj represents the HQ of habitat type j of grid x, and the value is in the range of 0-1; a high value means high quality. H j refers to the suitability (0-1) of habitat j. D xj z is the threat level of grid x, also considered habitat degradation. K is the half-saturation value, and z is a normalized constant. In our study, habitat suitability, habitat sensitivity to threats, weights and effective distances of degradation sources were all obtained from the literature [42,50,60,64,65] (appendix 1, tables 1 and 2, which can be found in the supplementaty information at is available online at stacks.iop. org/ERL/14/125006/mmedia).
Despite the abovementioned advantages in HQ evaluation, the biggest disadvantage of InVEST is that the threats outside of the study area boundary are always ignored. To minimize this disadvantage, we considered threats both within and outside the study Biodiversity hotspot analysis A broad definition of 'biodiversity hotspots'refers to any area or region with exceptionally high biodiversity at the ecosystem, species and genetic levels [51]. Spatial autocorrelation analysis can be used to identify biodiversity hotspots [43,66]. The * G i index is the most commonly used coefficient for spatial autocorrelation analysis. It is based on partial spatial autocorrelation using a distance weighted matrix, which can detect aggregates of high-value areas and low-value areas. The aggregates of high-value areas is called hotspot. The * G i index was proposed by [67], and it is expressed by the following equation: In this equation, x j represents the attribute value of the spatial unit in a partial spatial area. W ij is the distance weight between unit i and j; W i is the sum of all the distance weights; x is the average of all the attribute values in the study area.
i the value of the unit is higher than that of the neighbor unit; (with 90% or 95% confidence level), it is a high-value area, namely, a hotspot area [68]. When * < - (with 90% or 95% confidence level), it is a low-value area, namely, cold spot area.
In our study, we used ArcGIS spatial autocorrelation tools to conduct the hotspot analysis. For the habitat quality and NPP indicator, we used 90% and 95% confidence level as the high level hotspot and less high level hotspot. For the entire biodiversity hotspot analysis, we firstly kicked out the area where habitat connectivity is under 0.5, than we combined the hotspot area of NPP and habitat quality as the biodiversity hotspot. We classified the whole hotspot area into 4 classes. Namely, the biodiversity security pattern was classified into four classes. That is, when NPP and habitat quality value is both in high levels, we classified this area into first class biodiversity hotspot. When either NPP or habitat quality value is in high level, but the other is in less high level, we classified this area into second class biodiversity hotspot. When either NPP or habitat quality value is in high level, but the other is not hotspot, we classified this area into third class biodiversity hotspot. And the other hotspot area was classified into fourth class biodiversity hotspot area.

Urban growth simulation
Urban spatial expansion depends on five explanatory variables: accessibility, neighborhood, natural factors, planning, and socioeconomic variables ( [39,69], appendix 2: table 1 in the supplementary information). The construction of the urban growth simulation model with biodiversity constraints included the following steps ( figure 2): (1) quantify the influencing factor of urban growth and build the spatial logistic model. Obtain the potential spatial distribution of urban growth dependent on the size of urban growth from 2000 to 2010.
(2) Consider NPP, habitat quality, habitat connectivity, and their combinations as key biodiversity constraint indicators according to regional characteristics. (3) Take the minimum biodiversity loss as the total objective and simulate urban growth with multiple biodiversity constraint scenarios based on both the biodiversity loss and the potential spatial distribution laws of urban growth. For the details of the urban growth simulation, please see appendix 2 in the supplementary information.
In addition to the comparison of biodiversity losses among different scenarios, landscape pattern indexes were introduced to describe the spatial shape of urban growth. We used the contagion, connectance and aggregation indices to measure the condition of new urban spatial growth with or without biodiversity constraints. Then we used the connectance and aggregation indices to compare urban growth features in four directions (northeast, northwest, southwest and southeast) under different biodiversity constraints. The landscape pattern indexes were measured in FRAGSTATS 4.

Data sources
Accessibility to different centers is an important explanatory variable in the urban growth simulation (appendix 2 in the supplementary information). This article selects Tiananmen as the city center. According to the basic results of the second economic census yearbook in Beijing, Beijing Central Business District (CBD) and Financial Street were selected as the CBD area; Chaoyang CBD, Yizhuang, Tongzhou, and Jiuxianqiao were used as employment center; and Yizhuang, Shijingshan, Shangdi, Zhongguancun, Jiuxianqiao were selected as industrial center. To accurately portray the accessibility of the whole city, this paper will also fully considered the impact of urban subways on city commuting and used the vector raster incorporation method to obtain the comprehensive regional transport accessibility of the CBD, employment and industrial centers, subway stations, and high-speed tracks. The spatial data for the urban roads selected in this paper were derived from the transportation map of the 'Comprehensive planning of Beijing   The land use data were from the land use investigation data from 2000 to 2010. The socioeconomic data were derived from the Second China Economic Census, Beijing. The NPP data were obtained from the project of 'China's national ecosystem remote sensing survey of the year 2000-2010' [55]. The digital elevation model data were from the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences. All spatial data are presented in grid format, with a resolution of 100 m.

Results
Biodiversity security pattern analysis Features of biodiversity security levels From a biodiversity constraint indicator calculation, we obtain the spatial distribution of the biodiversity constraint indicator (figure 3). Using these maps and according to hotspot analysis, we obtain the spatial distribution of different biodiversity security level areas in the study area ( figure 4). Generally, the western area is the most important biodiversity security area. The area of first-class biodiversity security is 2871 Ha, accounting for 6.03% of the whole biodiversity security area. The majority of the firstclass area was distributed in the northwest area of Beijing. The first-class biodiversity security area is the core area of Beijing Xishan National Forest Park and Beijing Lufengshan National Forest Park; this area has a high habitat quality and a high NPP level and is far from most kinds of human activities (threats). This area consists of a nature forest and is well protected. The second-class area is 4084 Ha, accounting for 9% of the security area, most of which is around the core area of the national forest park, with a high habitat quality and a lower NPP level. However, this area is closer to the human activity area than the first-class area. The majority of the second-class area consists of nature forest but includes some shrub area as well. The third-class area has the largest area (27 382 Ha) of all the biodiversity security areas, accounting for 57.54% of the whole security area. This area has either high habitat quality or a high NPP level. The majority of the third-class biodiversity security area is concentrated in the western area of Beijing; meanwhile, its secondary largest distribution area is around the 6th ring road. This area consists of forest, grassland and tree nursery areas. The area of the fourth biodiversity security class is 13 250 Ha, accounting for 27.84% of the whole biodiversity security area. This area is more decentralized. However, this area has 3 major distribution areas. One of the major distribution areas is the transitional area from the national forest park to no protection areas; although it has a relativity higher NPP, it has a relativity higher threat from human activities. Another of the major distribution areas is the Yongdinghe River area, which consists of the water area and mud. Although the NPP level of the Yongdinghe River is low, it is an important area for aquatic beings. The last major distribution area is in the southeastern area of Beijing city, which is the crop cultivated area, with a high NPP level.

Habitat degradation analysis
According to the theory of InVEST, habitat degradation means the ecological land's impact degree by human activities. From the habitat degradation map of 2000, the degradation has a radial pattern around the core urban area (figure 5). The degradation degree decreased gradually as the distance from the city center increases. The degradation pattern is well in line with the urban construction and ring road distribution patterns. Specifically, we classify the degradation into 10 types according to Jenks natural breaks classification method [70]. The degradation degree of the strong impact area is 8.71%-12.98%, accounting for 9.89% of the whole ecological area (table 1). Most of the strong impact area is within the range of or around the edge of the 5th ring road. This area is the existing and intensive urban construction area, which has a very strong threat on ecological land. Moreover, the ecological land is sensitive to the urban activities. The degradation degree of the middle impact area is 4.07%-8.71%, accounting for 30.95% of the whole ecological area. The distribution pattern of the middle impact area is relatively disperse; however, a large amount of the middle impact area is distributed from the 5th ring road to an approximately 6 km distance to the 5th ring road. The impact degree of the light impact area is 0.00%-4.07%, accounting for 59.16% of the ecological area. Most of the light impact area is distributed around the 6th ring road.
According to the statistic of the average degradation degree of different ecological land use types, the degradation of the natural land use type is higher than the semi-natural land use type (table 2). The water area has the highest impact degree by threats, with an average degradation of 8.42%, almost falling into the strong degradation range. The reason for the strong degradation is that urban construction is more intensive around the water area, and high level roads tend to be built along the river. The lowest degradation land use type is cultivated land, with an average impact degree of 2.99%, which is in the light degradation range. Less construction in the cultivated land area is the most important reason for the low degradation. Another low degradation land use type is unused land. Most of them are distributed on the west edge of the study area, far from the urban construction area, thus less impacted. The other ecological land use types (forest land, orchard land and grassland) are in the middle degradation range.

Urban growth impact on biodiversity
Biodiversity loss under urban growth Generally, under different biodiversity constraint scenarios, all biodiversity losses increased gradually over time ( figure 6(a)). However, under different scenarios and urban growth phases, biodiversity loss characteristics depended on the value distribution character of the biodiversity indicators. The loss of NPP showed relatively smooth increasing trends during the study period. The reason is that the low value (under 0.5) NPP areas were relatively large and value variations were relatively small. Thus, under the NPP constraint scenario, urban growth occupies the low NPP value area gradually. The loss trend of habitat quality was relatively smooth before the 4th iteration but suddenly had a high loss in the last phase. According to the statistic of the habitat quality loss value, before the 4th iteration, the value of the mostly loss grid is under 0.5. At the last phase, urban growth searches for the high habitat quality value area. Regarding habitat connection, most of the low connection areas were within the 5th ring road range; these areas were the first to be occupied by urban growth. Then urban growth searches for the high connection areas; the connection value of these growth grids are almost as high as 1. Therefore, the connection loss after the 2nd iteration is relatively much higher; the connection loss percentage after the 2nd iteration is approximately 4%. As to the integrated biodiversity, its grid value is much affected by the habitat connection grid. Therefore, its loss trend is similar to the habitat connection factor.
Accumulative biodiversity losses of urban growth with biodiversity constraints is approximately 3% less than actual urban growth without biodiversity constraints ( figure 6(b)). The accumulative biodiversity losses of NPP are lowest either with or without the biodiversity constraint. However, when habitat connection is taken as the main constraint factor, much more biodiversity will be lost.
We made a correlation analysis between the habitat degradation degree and the biodiversity loss of the constraint factors. We found that the degradation's relevance (R 2 ) with habitat quality, NPP and the habitat connection area are 0.50, −0.12 and −0.09, respectively. The result means that when a habitat is affected by more human threats, the possibility is greater that it would be turned into an urban area. However, degradation has no significant correlation with the loss of NPP and habitat connection.
Urban growth occupation in biodiversity hotspot areas Urban growth from 2000 to 2010 occupied 14.80% of the biodiversity hotspot area (table 3), which is generally much more than it was with each biodiversity constraint. In the actual urban growth situation, 16.84% of the third-class area and 17.09% of the fourth-class area had been occupied. When referring to the biodiversity constraint situation, urban growth with the habitat quality constraint occupied the minimum (6.40%) hotspot area, with few firstand second-class hotspot areas, 6.86% of the third-class area and 8.54% of the fourth-class area. The hotspot area occupation is similar under the other 3 constraint factors. In that case, the third-class area is largely occupied.

Urban landscape patterns under different growth situations
Compared to the landscape patterns of the new urban growth area, the actual urban growth area has more sprawl, more aggregation and less connectivity than urban growth under the biodiversity constraint. As showed by the landscape change pattern under different urban growth situations (table 4), the urban sprawl value of the actual urban growth patches (49.77%) is much higher than urban growth under the biodiversity constraint (approximately 25%). Meanwhile, the connection degree of the new growth patches in actual urban growth (0.33) is lower than urban growth under the biodiversity constraint (approximately 0.5). Finally, the aggregation index value in actual urban growth is 63.52%, much higher than in the biodiversity constraint situation (below 57%). Urban growth under the biodiversity constraint is relatively split because the lower biodiversity area is relatively disperse; urban growth priorities occupy these areas under the biodiversity constraint.
We also analyze the connection and aggregation features of new urban patches at different directions under the biodiversity constraint. We observe that the connection degree of new urban patches is much higher in the southeast, with a connection value of approximately 1 ( figure 7). The second highest connection area is in the northwest. In the northeast and the southwest, the connection values are relatively lower. The aggregation pattern in the northwest is much lower than it is in the other 3 directions. Moreover, the aggregation pattern in the southeast is much higher than it is in the other 3 directions. Particularly, in each direction, the aggregation pattern under the NPP constraint is higher than it is in the other biodiversity constraint indicators.

Discussion
Biodiversity constraint indicator for urban spatial optimization Biodiversity conditions reflect ecosystem health [71]. Biodiversity has a strong synergistic relationship with  ES [40]. As urban growth optimization needs to maximize ES [23], building a biodiversity security pattern could be a priority of urban development. In our study, we establish a biodiversity constraint indicator using urban growth constraint factors, including NPP, habitat connectivity and habitat quality, which makes efficient control of biodiversity loss as urban growth. Combining these 3 indices, we generate a biodiversity security pattern in the study area as well. These 3 indices are crucial for keeping biodiversity and the ecosystem stable and integrated, as these factors represent high-level ES, habitat connectivity and habitat degradation at the same time. As many studies tend to choose different kinds of ecological factors in urban optimization, which makes regional application and comparisons more difficult [16,[18][19][20]39], establishing the scientific biodiversity constraint indicator as a basic constraint could be a good method of maintaining urban ecological security.
Studying an interaction means simultaneously considering urban growth and biophysical factors and providing important feedback mechanisms [3], which can also help generate a theory of interactions and feedback mechanisms between human decisions and ecological resilience in urbanizing ecosystems [72]. Most currently urban planning methods simply treated ecological factors as special factors rather than as a comprehensive factor [11]. As city growth and the urbanization process are strong linked with biophysical and ecological processes [1], optimal urban expansion simulation should combine ecological factors with the urban growth process. In our study, urban growth simulation firstly considered the future potential spatial distribution area based on five major aspects: accessibility, neighborhood, natural factors, planning, and socioeconomic variables. Then simulation took minimum biodiversity loss in the iteration of urban growth. Therefore, ecological factors are interacted with other urban growth factors. The InVEST model was used to obtain habitat quality in our study. In this model, habitat quality depended on the ecological land and human stress, which included human and ecological interactions. Moreover, we consider the stress both within and outside the study area, namely, all of Beijing city, which makes the interaction more accurate.
Future studies could make some improvements to the application of the biodiversity constraint indicator. For urban growth simulation, we can use the biodiversity constraint indicator to obtain the biodiversity hotspot first and make it the restricted area for urban growth. Therefore, we can maintain the biodiversity security area not being occupied by any construction. We also need to improve the calculation of the habitat connection because we hypothesized that all kinds of habitat (e.g. forest and agriculture land) had the same capacity for connectivity. Therefore, habitat connectivity in most areas is relativity high, which aligned the integrated biodiversity constraint indicator with the habitat connectivity value. A more accurate integrated biodiversity constraint indicator can be better for urban growth optimization than separate biodiversity constraint indicators.
Simulation results' implications for future urban development The distance from the city center will have an important impact on urban expansion. As the distance from the city center increases, the possibility of land development will be significantly reduced, which will directly lead to the urban sprawl structure of blind expansion in Beijing [73]. Existing findings also indicate that the distance from the city center is a primary factor that affects urban expansion in other cities [74][75][76]. To prevent urban sprawl [77], the coordinated development of Beijing-Tianjin-Hebei and new city subsidiary-center of Beijing have been successively implemented in recent years. As the new center of Beijing, Tongzhou will become an important area where high-end resources will rapidly gather in the process of accelerating the reconstruction of the capital, which is bound to have a new impact on urban expansion in Beijing. The layout of the city center, which resulted from administrative power, will greatly affect the speed and quality of urban expansion [78]. Such short-term breakthroughs in the traditional inertia effect and development pattern will lead to rapid changes in urban morphology, especially urban sprawl caused by rapid expansion [79]. New city subsidiary-center construction cannot form an agglomeration effect or be finished as planned because an accurate prediction of urban development needs and a scientific demonstration of the urban center location are lacking [80,81].
The existing urban growth control policy guides the orderly growth of urban space and ecological environmental protection through an accurate prediction of urban growth. However, prediction mistakes will invalidate subsequent control policies, which will lead to new urban growth and ecological problems [39]. In this paper, a biodiversity priority-oriented urban growth simulation method is proposed. This proposal is an optimization plan for urban growth under the guidance of 'ecological forced'. Under the guidance of ecological optimization, urban growth in the research results reflects decentralized and multicenter spatial development characteristics. Furthermore, this type of urban growth not only provides a new model for breaking the inertia of urban sprawl but also proposes 'ecological security' as an applicable regulatory tool for urban planning and space governance reform.

Conclusions
Facing the limitations of selecting more scientific ecological indicators for urban expansion optimization, we establish a basic biodiversity constraint indicator that includes NPP, habitat connectivity and habitat quality as the biodiversity constraints to integrate urban growth simulation. Based on these indices, we found the biodiversity hotspot in the study area and efficiently controlled biodiversity loss for urban growth. We also found that the habitat degradation pattern aligned with the urban construction and the ring road distribution patterns. Threats from the urban area and the transportation system are the main reason for habitat degradation in city areas. Moreover, degradation of the natural habitat type is much higher than that of the semi-natural habitat type. Finally, under the guidance of the biodiversity constraint, urban growth reflects decentralized and multi-center spatial development characteristics, which provide the potential to break the inertia of urban sprawl. This kind of urban growth under a biodiversity security pattern can be an applicable regulatory tool for urban planning and space governance reform.