When Shape Matters: Using a Simple Mathematical Model to Estimate Critical Area Sizes in Conservation


 In the analysis of anthropogenic impact on the environment arises the question of whether the shapes of preserved habitat fragments play an important role in the conservation of wild species. In this work we use a very simple mathematical model based on a reaction-diffusion equation to analyze the effects of geometric shape and area on the permanence of populations in habitat fragments. Our results indicate that a dimensionless quantity calculated from a combination of biological variables is the main component that determines if the species survives in the preserved fragment and whether its geometric shape is important. We provide a methodology to calculate critical area sizes for which population size is most affected by fragment shape. The calculation is based on four quantitative variables: maximum per capita reproduction rate, per capita mortality rate outside the fragment, carrying capacity in the conserved environment and mobility in the disturbed environment. The methodology is illustrated by a preliminary study, in which the model is used to estimate threshold area sizes for habitat fragments for the threatened species Sapajus xanthosternos .


INTRODUCTION
Landscape ecology studies traditionally use landscape pattern indices, a.k.a landscape metrics, to predict ecological responses, they are mostly associated to patch size, shape and habitat amount and aggregation 1 . Detailed ecological research on patch-scale alterations on biodiversity were never as relevant, considering that half of the forests of the planet are less than 500 m from the forest edge, and mostly are patches smaller than 10 ha 2 . Moreover, habitat fragmentation, isolation and creation of edge environments initiate long-term responses of organisms and ecosystems processes that percolate through the landscape 2 .
The edge of forests may be a barrier to animal movements 3,4 , thus, the area and shape of forest patches are appropriate metrics to assess the effects of spatiotemporal changes in the landscape configuration over biodiversity measures, such as, species richness, community structure and organisms' abundance 5,6 . These alterations in landscape patterns are also used as predictors of ecological processes, namely the probability of population extinctions and migrations 7 .
Recently, a vivid debate on the importance of considering several spatial scales to assess the changes in biodiversity has arisen and most studies indicate that scale selection is species-sensitive 8 . It is well-known that anthropogenic and natural landscape alterations disproportionately impact forest-core species, which are almost four times more prone to extinction than edge-tolerant and habitat generalist species 9 . Accordingly, it has been shown that forest-core species are vulnerable to hunting and predation when moving through non-forest habitats, and that abundance of forest-core animals is consistently larger at about 400 m away from the edge, confirming that edge effects operate at a small spatial scale 9 .
The impacts of edge effects on populations of forest-dependent species are influenced by the ratio of forest core in relation to forest edge. Thus, edge effects are strongly related to the shape of forest fragments, but they are less intensive in the core of larger forests fragments 10,11 . Forest-core species are more strongly affected in small forest fragments with convoluted shapes 12,11 .
The population abundance of these forest species tend to decline under edge effects and ultimately increase the risk of local extinction 9 .
Traditional approaches in landscape ecology use diversity measures, such as species richness and equitability, to explore how changes in habitat structure and landscape configuration may lead to species loss, alterations in the animal community composition, and biodiversity erosion 13 . Fewer landscape ecology studies address the influence of those landscape modifications on population dynamics. Here we use population dynamics models (partial differential equations) to understand the population vulnerabilities associated to changes in the spatial configuration of forest patches. This approach allows us to study the link between the spatial configuration of the habitat and the population dynamics 14 .
The paper is subdivided in three other essential sections besides this introduction. In section 2 we present the model, the biological hypotheses and parameters, and quickly review some metrics of fragment shape. In section 3 we deduce general ecological results from the model and illustrate its application in a particular preliminary study directed to a particular species (Sapajus xanthosternos). Finally, in section 4 we discuss some of the general and particular results, delineating also directions for future research.

The model
The main focus of this paper is to deduce conditions for which the shape of a habitat fragment does have a significant impact on the chances of conservation of a certain species that occupies it. Since this question is quite general, we are immediately caught in a tension, trying to avoid two opposing dangers: oversimplification of biological traits and narrowness of results due to an excessive number of hypotheses.
For the particular model used in this work, each species may be represented by a combination of five basic biological parameters: the time it would take a population to be doubled under the most favorable conditions ( in time units), the carrying capacity in an undisturbed/preserved environment ( in number individuals per area unit), the time that it would be necessary for the population to be halved in a disturbed/hostile environment ( in time units), a dispersion coefficient ( units of area −1 per unit of time) and a non-dimensional coefficient of mobility ( ) that identifies if individuals tend to disperse faster or slower in the disturbed environment when compared to its movement in the preserved areas.
Of course, under many aspects, this is an oversimplification of any biological species, but here we argue for this approach, underlying the following points: 1. If our objective is to gain insight, from a broad perspective, on the effects of shape and area size on the survival of species, we cannot rely on detailed dynamics which are species-dependent.
2. This generality allows for an easier application in specific cases, requiring just a few biological parameters to adequately apply the methodology outlined by the model. Detailed models often require the estimation of a much higher number of parameters, many of which may demand intricate experiments to be successfully determined.
3. As long as the model and methodology is not considered as an absolute tool and is used only as a guideline to provide rough estimates of threshold area sizes, there should be no harm in approximating the complex population dynamics of the real populations for the simplified one presented here.
It is worth to mention that the model is effective in illustrating the importance of the relations between certain scales that are connected with mobility, reproduction rates and area of the fragment. In this sense, the model should provide some insight from the perspective of ecological theory and conservation, going beyond the simple applications of statistical methods of analysis, which, of course, have their own importance and scope.

The biological hypotheses
In the process of mathematical modeling it is important to clarify the biological hypotheses assumed so that we know exactly what is included and what is not in the dynamics presented by the model. The main assumptions considered in our model are enumerated below: 1. The model describes the dynamics of one particular species at a time. No interspecific relations are modeled.
2. The region to be analyzed can be divided into two clearly different landscapes. One represents a preserved and favorable environment to the species while the other corresponds to a hostile ambient, where the population could be sustained only for an definite amount of time (extinction is unavoidable in the hostile ambient).
3. Within the preserved environment the species grows logistically up to a carrying capacity. 4. In the hostile ambient the population decays exponentially.
5. Individuals are forced to move around in the environment in both types of landscapes, either due to overpopulation forces, foraging or other species-dependent factors. The particular details are not taken into account in the model, instead, those forces are represented by two parameters, one for movement in the preserved region and other for the hostile one.
The first hypothesis simply means that we are not explicitly describing complex ecological relations between species. These relations, of course, are very important in the real biological system and may be the decisive factor between the survival or the extinction of a particular population. In our approach, those relations are "embedded" in the form of the division of the region in a preserved or hostile environment. In the preserved region, all the ecological relations necessary for the survival of the population are present while in the hostile region those conditions are absent. In this way, the analysis can be focused in just one species.
The second hypotheses is related to a very frequent occurrence in the study of anthropogenic impact on wild populations. In many instances is is possible to observe a "natural" landscape, where wild populations reproduce and a "disturbed" one, that has been modified to serve another purpose, for instance, raising of cattle or agriculture 15 . Although the biological hypotheses do not explicitly account for a continuous transition between those types of environment, in the sense that each location is either classified as "preserved" or "disturbed", the mathematical formulation does provide a certain smoothness in the transition from one scenario to the other. In this type of model, populations that live closer to the transition frontier tend to reproduce, on average, slower than those in more "interior" areas.
Referring to the third hypothesis, although some species tend to conform to the specific form of Logistic growth 16 , in our model, this particular form of dynamics is used just to represent the factual observation that every population has its growth limited to the natural resources available. In this case, those are supposed to be directly proportional to the size of the preserved area. Similar models with analogous limited-growth functions (a generalized Logistic model was used 17 ) were tested and very close results were obtained so we have some evidence that our analysis in independent of the particular choice of a Logistic growth, but rather that the fundamental fact is the limited carrying capacity of the environment.
The fourth hypothesis accounts for mortality in the hostile environment. It can be proved 18 that if each individual has a fixed probability of dying per time unit then the dynamics of the whole population can be approximated by an exponential decay. The important factor here is that the species is tending to extinction in the hostile environment.
The last hypothesis involves the dispersion of the individuals. Animal movement is an extremely complex topic and many factors do contribute to the observed movement pattern at the population level. When enough data is available, precise mathematical models can be fitted to describe complex inter and intra-specific interactions. In particular, Moorcroft et al. 19 have successfully used a reaction-diffusion equation to model wolf pack territory formation in Yellowstone Park, being able to closely fit the model to empirical data on population distributions. Due to its generality, the model we propose is much less detailed and we do not expect such model to precisely describe the microscopic individual movement dynamics. Instead, we rely on a macroscopic parameter, the diffusion (or dispersion) coefficient, to represent a general tendency of dispersion, whatever is the particular biological mechanics behind it. Particular population estimates in well-known regions may then be used to fit the dispersion coefficient, avoiding the need of intricate empiric research on the movement pattern of the modeled species.(see Section for further details).
Finally, there is the observation that individuals may move at different speeds depending if they are in a preserved or hostile area. While the precise estimation of such difference may be subject to technical difficulties, we find this assumption necessary, given the generality of the model, which could be applied to many different species and terrains. In Section 3.1.2 we show that parameter , which is related to this difference in movement, is a parameter that does not have a very strong impact on the equilibrium population. Since the results are fairly robust to this factor, the model may still be applicable even if it is not possible to estimate precisely the difference in mobility between hostile and preserved areas.

The variables and the simulation space
The dynamics proposed occurs in a three-dimensional space (ℝ 3 ), having two Cartesian coordinates, ( , ) ∈ ℝ , for the spatial distribution and another one for time, ∈ ℝ. The population is then described by a density function ( , , ), so if Ω ⊂ ℝ 2 represents the fragment of preserved habitat while − Ω stands for the hostile environment. Figure 1 illustrates an example of this type of separation.

The reaction-diffusion equation
Reaction-diffusion models are applied to the study of a wide variety of natural systems. Pattern formation in animal coats 20 , spatial distribution of slime molds and formation of galaxies 21 , ecological invasion by alien species 22 , chemical signalization in insects 23  To model the population dynamics in , we propose the following equations: where and are the dispersion coefficients in the preserved habitat and the hostile environment, respectively, is the per capita reproduction rate in the habitat, is the per capita mortality rate in the hostile environment and is the carrying capacity in the preserved habitat.
The dynamics is separated in two regions, Ω and − Ω. The first term in both equations represents the mobility of the population in each environment, so the ratio = ∕ indicates where the movement is faster, if > 1 individuals move faster in the hostile environment while < 1 means the opposite.
The use of a diffusion equation to incorporate the characteristics of complex movement patterns may be criticized as an oversimplification. While it is true that this type of equation may be deduced in many different ways 26 , a common misconception is to think that it can only represent Brownian motion, that is, the movement displayed by particles randomly colliding and moving in every direction. One way to interpret differently this equation is to suppose that, in microscopic scale, individuals carefully avoid overcrowded regions, moving in the direction of less populated areas in a velocity that is proportional to the stimulus given by the Weber-Fechner law 27 . Under such assumptions, the resulting macroscopic dynamics is exactly a diffusion equation [1] .
This can be seen by simply writing the term corresponding to the flux of individuals, The biological interpretation is clear now: individuals move on the opposite direction of the gradient, but with a a speed/probability that is proportional to the ratio of the concentration difference perceived and the concentration of individuals. This ratio incorporates the effect described by the Weber-Fechner law, that the perception of changes in the intensity of stimuli is relative to the total present stimulation. Just to give a simple concrete example of this psychophysical law: a match lighted in total darkness provides a much higher stimulus than one lighted in plain daylight. So the model can be interpreted as implying that the population is trying to avoid overcrowded areas, moving in the opposite direction of increasing population density. Although the Weber-Fechner is not an "absolute" natural law and has been subject to criticism and also that alternatives formulations have been proposed for modeling stimuli perception 28,29 , recent results seen to find new evidence supporting it 30 . In any case, here it serves the purpose of illustrating how the diffusion equation may be re-interpreted as portraying a more complex behavior than just Brownian motion.
As stated in section 2. 1. A rectangular domain facilitates the implementation of numerical methods to solve the differential equation.
2. The domain is chosen in a way so that it included that fragment of interest (Ω) and that the boundary conditions do not have a significant impact on the results.
The boundary conditions for the mathematical problems may be homogeneous Dirichilet conditions: where stands for the boundary of . Another possibility is to use the Robin homogeneous condition is given by: where ∕ is the normal derivative (pointing outwards region ) and > 0 is a constant that defines the intensity of the flux.
Homogeneous Dirichilet conditions imply that individuals die at the boundaries while Robin homogeneous conditions represents a migration that is proportional to the population density at boundaries We will not stress too much the role of the boundary conditions because the focus of the model is to analyze the effect of the geometric shape and area of Ω. The domain , in this case is just a convenient choice that allows for an easier numerical simulation. We tested our results using both conditions and a variety of values for and the results were very similar, independently of the boundary conditions used..
One last technical observation on the model is that there is a discontinuity in the diffusion coefficient. This is not a relevant modeling problem for our particular case since, in this work, the solution of the differential equation is approximated using a discretized version of the domain, where discontinuities are not an issue.

Numerical methods
To simulate the model we first divided the domain into the two regions Ω and − Ω and, observing that and are constant, obtained the two partial differential equations: where Δ stands for the two-dimensional Laplacian operator.
To approximate the solution of those equations, we adopted a finite-difference scheme (centered, second order) for the spatial Laplacian operator (Δ) and Runge-Kutta methods of combined orders 2 and 3 to solve the resulting system of ordinary differential equations. The application of these methods is considered routine in numerical analysis and can be consulted, for example, in Burden and Faires 31 .
To implement such numerical schemes we used two softwares: Matlab R2016b (registered for academic use under license number 1115837) and the open source Scilab (version 6.1.0). Two independent teams worked on the codes, each using a different software for comparison of results and avoidance of coding errors. For the Runge-Kutta methods we used the built-in functions provided by both softwares: ode23 for Matlab and ode (with method option "rk") for Scilab.

Estimates for , and
The model uses five mathematical parameters to simulate the dynamics: , , , and . To estimate those parameters, we need some biological estimates, as mentioned in Section 2.1: 1. : supposing the conditions are very favorable to the population (resource abundance, low competition), is the time the population takes to double its numbers.

2.
: supposing the whole region was transformed into a hostile environment, is how long would it take for the population to be halved.
3. : if we denote by the average daily displacement of individuals in the preserved habitat and in the hostile environment, then may be interpreted as an estimate of ∕ . Here the time scale of one day was mentioned, but it can be adapted according to the biological convenience.
The exact method for providing those estimates are left to the expert biologists for each case. From these three biological parameters and some population estimates in known areas and a least squares fit it is possible to estimate the remaining parameters ( and ). In Section 3.2 we provide a case study as a guideline.
Parameter : From the estimate of , and taking into account that, for populations that are small relative to the carrying capacity, the Logistic growth can be approximated by a Malthusian one, we obtain as: Parameter : Since the population decays exponentially on the hostile environment, the instantaneous rate of decay can be directly obtained from the estimate :

Measures of compactness
One of the central questions approached in this paper is to evaluate how important fragment shape is to the conservation of populations. To conduct a mathematical analysis we need, then, some measure of "compactness" that can be used to represent the characteristics of a certain shape, so that its impact on population can be quantified. A discussion of the roles of such measures and a review of some particular formulas can be consulted in Li et al. 32 . Other authors 33,34 present also some shape measures commonly used in landscape ecology analysis. Below, we briefly review some of the measures and choose which suits best the scope of this work.
If Ω is the two-dimensional set that represents the fragment we denote by (Ω) as the perimeter (we suppose that the boundary of Ω can be described as a smooth function) and (Ω) its area.
Mcgarigal et. al 34 [p.104] present two of the most common measures of compactness, PARA: and a Shape Index (SHAPE): Another possible measure is the IPQ 35 , given by: The measure given by PARA is not convenient because it is scale-dependent, so it was not adopted in this work. SHAPE is not dependent on scales and it is directly related to the IPQ.
In the IPQ measure a circle has the maximum measure of compactness value of 1. The intuition is simple, as the perimeter increases for the same area, the measure decreases, resulting in a value between 0 and 1. The use of the square of the perimeter avoids the impact of scaling factors (units of length), resulting in a reliable measure. To make the results more intuitive, we use the acronym GE to stand for "geometric efficiency" as the measure of compactness.
Here we propose a slight modification of the IPQ, using for the measure of compactness the following formula: The change is only a re-scaling of IPQ but, in our simulations, the discretization of the domain causes the regions to be transformed into a union of small squares, which also affects the length of its perimeter and, to a minor degree, its area. Due to this slight deformation of regions, we find it useful to use the square as the reference figure instead of the circle. In this measure, the square has a compactness measure of value 1.
SHAPE is directly related to GE: so any results obtained in terms of GE can be easily converted to SHAPE and vice-versa.

RESULTS
This section divided in two parts. The first is dedicated to general results and the second to a case study of the species Sapajus xanthosternos.

General results
For a clearer presentation, we separated the general results into three further sections. First, the combination of parameters that are important to the simulations are identified, then a sensitivity analysis is performed and finally the impact of geometric shape on species population is discussed.

Identifying the key parameters
Before we addressed the question relative to the impact of geometric shape of a fragment on populations, it was important to analyze the role of the parameters when shape was fixed. We chose a fragment as a circle and used six parameters to simulate This re-scaling provides an estimate of the border effects on the fragment. For instance, if * = 0.8, this implies that the influence of the outside hostile environment causes a decrease of 20% of the maximum potential population in the fragment. By simulating the model it was possible to obtain the re-scaled population as a function: * ( , , , , ) and we used this relation to study how the parameters affected the survival of the species. Allowing all parameters to vary in each simulation, the trend was not clear * when we looked individually at parameters , , , or , while parameter had no impact at all in * . Instead, we found that * is dependent on the non-dimensional groupings: and the non-dimensional parameter .
To illustrate such dependency, we chose three different values of (1∕2, 1 and 2) and for each value of , the following distributions for the other parameters was chosen: where  ( , ), , ∈ ℝ stand for a uniform distribution between and . Note that once and are established, and are automatically determined by the relations = and = .
These intervals were conveniently chosen to illustrate the relations and similar results are obtained with other intervals that comprise a similar spectrum for . It is worth to stress that since , and are non-dimensional, they are not affected by the choice of units used to describe population density. A total of 250 random combination of parameters was chosen for each fixed value of . In Figure 2 we present illustrations of how * behaves in relation to some of the parameters, this illustrates the important role of and . The results were robust to the choice of both initial and boundary conditions.
Another way to analyze the impact of the parameter is to measure the correlation between the parameters and * . In Table   1 we present the Spearman coefficient 36 of correlation for the parameters and the variable * . The results confirm that is the most important factor in the determination of * .

TABLE 1
Spearman correlation coefficient for the simuations of this section. , has the largest value, illustrating its importance to the survival of the species. Since parameters , and are directly related to , there are also a significant correlation.

Spearman correlation coefficient for
wherē andΩ are transformations of the domain and the fragment region Ω to ( , ) mapping. For the particular case of the simulations in this section, Ω was taken as a circle of area having a radius of = √ ∕ , so the transformed regionΩ is a circle of radius̄ = ∕ √ ∕ = √ ∕( ). Re-arranging the terms, we obtain: Equations 16 and 17 show that the behavior of the non-dimensional model (which is independent of scaling) depends only on the non-dimensional groupings , and (which shall be called parameters henceforth).
The results of both the simulations and the dimensional analysis show clearly that parameters and are essential to the survival of the species. These results are biologically intuitive: 2 is proportional to both and , meaning larger areas and fast reproduction rates favor the survival of the species (hence the positive correlations between and * ). Also, 2 inversely proportional to , meaning that large mobility associated with a hostile surrounding environment impacts negatively the population. Parameter indicates how hostile is the environment outside the environment, being proportional to and inversely proportional to , with more hostile environments leading to lower populations. Parameter had a smaller influence on the final results, as can be seen in the curves in Figure 2 -d), for each fixed curve with fixed , was allowed to vary between 0.5 and 2, with little impact on the final re-scaled population * which was also confirmed by the analysis of the correlation coefficients in Table 1 .
To analyze the relative impact of parameters , and we performed a sensitivity analysis of model 16, presented in the section below.

Sensitivity analysis
In this section we consider * as a function of , and . The simulations are analogous to the ones in Section 3.1.1. We chose the shape of the fragment as a circle (with radius given by Equation 17), taking all parameters with the same distribution: To evaluate the sensitivity we used the slopes of the direct one-dimensional regressions ( * × , * × and * × , from the multilinear regression ( * = + + + ), the standardized regression coefficients and the Spearman coefficient. In Table 2 we display the results of the analysis. They clearly confirm the strong influence of parameter , with and playing a secondary role.

The impact of geometric shape
In the previous sections it was clearly established that and are key parameters that define species survival in the model. The simulations indicated that, for each fixed , it is possible to divide the region of variation of in three quite distinct regions. For very low values of , the population is unable to survive, with * → 0. For very high values, boundary effects have little impact on the dynamics, so * ≈ 1 with the population approaching the maximum value possible. Finally, there is an intermediate region where * makes a transition between those extreme values, and this is strongly related to the effect of geometric shape.
For very large or very small values of , geometric shape in unimportant, because either the area size is too big or too small in relation to the scales of reproduction and mobility (which is the measure that parameter incorporates). In the intermediate region the geometric shape plays crucial role in determining how strong is the boundary effect, which may lead either to extinction or preservation.
The central problem analyzed in this paper can now be addressed as a very precise question: for each point in the twodimensional parameter space ( , ), is it possible to determine if geometric shape affects significantly the population estimates by the model? To answer this question we designed procedures do create random shapes and also to deform a preexisting shape, affecting its measure of compactness (as defined in Section 2.4).
For each point in the ( , ) map we took an initial geometric shape for the region Ω, which we defined as 0 . This shape had a total area equal to 2 to correspond to a dimensional simulation of a figure of area and a initial GE, 0 = GE( 0 ) larger than 0.65, to allow for variation. By applying deformations on this initial shape we generated a set of 30 other shapes, 1 , … , 30 , with the same area but with compactness measures 1 , 2 , … , 30 in the set (0.05 ⋅ 0 , 0 ). In Appendix A we discuss in more detail the methods used to create deformations of the initial figure as to obtain shapes with smaller compactness values GE. This process is then repeated 10 times, so we obtain a collection of 300 shapes for each ( , ) point.
For this set of shapes we simulated the non-dimensional model for a time horizon of = 100 (which corresponds to the time horizon of = 100∕ in the dimensional model) calculating the final re-scaled population * = ( )∕ 2 for each shape, where ( ) denotes the total population inside after the simulation of the model. We obtained thus numbers between between 0 and 1 that correspond to the fraction of the maximum possible population achieved by that particular combination of parameters and geometric shape.
Given the series of points ( , * ) relating figures with different measures of compactness and the final re-scaled population, it is possible then to calculate regressions of * by that estimate the impact of the measure of compactness in the final rescaled population. Since both variables are non-dimensional and are scaled between 0 and 1, we adopted the slope of the linear regression to estimate the correlation, obtaining the impact of geometric shape as a function of and , ≡ ( , ).
In Figure 3 we present a scheme that illustrates the whole process described in this section.

FIGURE 3
Scheme for the simulations of geometric shape impact on population. Each initial shape is deformed, reducing its compactness, and the impact on population recorded.The distribution ( , ), of linear regression slopes, measures the relation between compactness and population levels. In Figure 4 we present examples of distributions of , * and a curve ( , ) for = 1. The distributions clearly illustrate that, when = 1, we have the maximum impact of the geometric shape for values of ∈ (5, 10).
In Figure 5 we present the map of slopes obtained from the simulations, ( , ). In the map is drawn also two curves ( ) = 0.1348 0,007490604 3 and ( ) = 0.0096 0.2494 . * is used to estimate area sizes that are greatly impacted by fragment shape. To estimate this curve, for each fixed value of we estimated the corresponding for which ( , ) was at it maximum value, next using a least squares fit, we adjusted a curve to fit the points. In this way * represents critical points where fragment shape has maximum impact and can be used to estimate areas sizes ( * ) for which fragment shape is very important. Finally, Please note that the relatively small effects produced by variations in are not included in these equations

A preliminary case study
In this section we present a case study for the species Sapajus xanthosternos. We must stress that this study has the main purpose of illustrating the application of the model and how to interpret its results. For more robust biological conclusions concerning the species, further empirical work in estimating the parameters should be conducted.
It is possible to estimate the maximum per capita growth rate if the maximum annual growth rate of the species is known.
For this particular species, it was estimated that, given the most favorable conditions, the population would have a maximum growth rate of 13% per year 37 , which leads to a correspondent logistic growth rate of = ln 1.13 ≈ 0.1222∕year.
Mortality rates in the hostile environment are harder to estimate, and we had to rely on an educated guess by the expert to arrive at a reference value. The information provided is that a group of 64 individuals could not last longer than a year, if exposed continuously to the hostile environment in this time interval. With such information and, considering an exponential decay in the hostile environment it is possible to obtain: = ln 64 ≈ 4.1588∕year.
Since Sapajus xanthosternos is an arboreal species and the hostile environment, in this particular case, is considered to be mostly deforested areas dedicated monoculture, individuals tend to move and disperse faster in the impacted than in the preserved areas. The meaning here is that even though the individuals might move faster in arboreal environment, they tend to move more often in deforested areas, always looking for forest cover, resulting in a higher average speed in the open areas. Considering this information and also by expert advice, we adopted then a value of = 2, meaning that, on average, individuals move two times faster outside the preserved areas.
To estimate the remaining parameters, and , we collected information on population estimates of populations of Sapajus xanthosternos 38 living in ecological sanctuaries near hostile environments. Since these population estimates are based on extrapolation from number of individuals sighted, we selected only those studies that were conducted in smaller regions (total area below 1000 ha) so that these estimates can be considered more precise. Ii is also worth to mention that 1000 ha is the largest forest area size in which Sapajus xanthosternos was recorded 39 . After this process we obtained 5 fragments with different area sizes, shapes and populations. Such data can be consulted in Table 3 .
To obtain images of the fragment shapes we used the software QGIS, using the CBERS4A downloader plugin, which provides access to the online database of CBERS4A satellite images. The images used correspond to the years in which the population estimates studies were conducted each individual fragment. These are displayed In Figure 6 . For each region, the model was simulated until the population reached an equilibrium. At the end of the simulation we calculated the total population inside the protected region, as in Figure 6 , obtaining thus a vector with five population estimates provided by the mathematical model.
We then used the least squares method to determine the optimal values of and that led to the best fit.   (5). Simulations of the model were performed in these domains and the results were compared with the field estimates for the populations (see Table 3 ).
In Table 3 we present the results from simulations using the least squares fit for and . The optimal value for was * ≈ 9.1746 × 10 −4 year −1 ⋅ km −2 and for was * ≈ 26.8033 × 10 −4 individuals ⋅ km −2 . This particular value of * was considered biologically plausible by the expert biologist and within his own estimates.
With the estimated value of * now it is possible to use equations 19 and 20 to establish to estimate the most sensible ( * ) and a relatively safe ( 0.4 ) area sizes for the species Sapajus xanthosternos. Results are presented in Table 4 .

DISCUSSION
The simple model developed and the simulations allowed us to identify clearly the most important parameters that determine how strongly populations are affected by fragment shape. The results show clearly that these parameters are associated with mobility of the species, reproduction rate, area of the fragment and mortality rate in the hostile environment. All these parameters are combined in the two non-dimensional groupings and .
The two-dimensional mapping of the behavior of the model in relation to parameters and allowed us to establish functional relationships (i.e. 19 and 20) for a * , area size for which the population is very sensible to geometric shape, and 0.4 , an estimate of a sufficient large are where geometric shape does not have much impact on population size. These relations shows is that the area threshold is proportional to the mobility (represented by ) and inversely proportional to (rate of reproduction), reflecting the intuitive ideas that more mobile species need larger areas and those who reproduce faster need smaller ones. Also, the fact that the parameter that represents mortality in the hostile environment ( or ) is present as an argument of a logarithm function reflects the fact, also shown in the sensitivity analysis, that its impact on the area threshold is smaller.
Another, somewhat obvious, biological result that can be derived from the equations is that species with smaller body sizes should also need smaller areas for their preservation. This is because smaller species tend to reproduce faster 40 , meaning that

Availability of data and material
All the relevant data used in the manuscript is available in the published scientific references. The software used in the simulations is licensed for academic user and is registered in the name of the first author.

Code availability
The codes used in the simulations are available upon request by electronic mail to the first author. We observe, however, that the authors do not make the commitment to provide an user-friendly code and/or explain the details of the implementations. The methods used in the codes are explicitly described in the manuscript.

Author contributions
Raul

A GEOMETRICAL METHODS OF SHAPE DEFORMATION
To generate figures with various values of GE, we adopted three methods. Method 1 is based on "stretching" the figure in one direction, method 2 is based on adding irregularities at the boundaries, increasing the perimeter without stretching the figure in any direction, finally, the method 3 is a combination of the first two. In figure A1 we present some examples of figures generated by the three methods.
Given a target area 0 , a polygon with vertices is generatec using polar coordinates, with angle = 2 +1 , , and ray , = 0, 1, ⋯ , − 1, such that the area of is 0 . For each region , two preserving area geometric transformations are applied The data used in this article was create according to the following steps: 1. Initially, for a given 0 , a region is created with vertices, following a uniform distribution ([ , ]) and following a uniform distribution ([ , ]) such that > 0.65; 2. Using the first transformation method on , ten regions are created with geometric efficiency 10 , = 1, 2, ⋯ , 10; 3. Using the second transformation method on , ten regions are created with geometric efficiency 0. 4. Using a combination of those two methods, ten regions are created like this: first,̃ is created, using the second transformation method, with geometric efficiencỹ = 0.2 + −1 12 and then is created, using the first transformation method, with geometric efficiency ̃ , < 1.
For a given 0 , steps 1-4 previously describe are repeated ten times so that, for such 0 , thirty regions are created.