Restoration Ecology: Two-Sex Dynamics and Cost Minimization

We model a spatially detailed, two-sex population dynamics, to study the cost of ecological restoration. We assume that cost is proportional to the number of individuals introduced into a large habitat. We treat dispersal as homogeneous diffusion in a one-dimensional reaction-diffusion system. The local population dynamics depends on sex ratio at birth, and allows mortality rates to differ between sexes. Furthermore, local density dependence induces a strong Allee effect, implying that the initial population must be sufficiently large to avert rapid extinction. We address three different initial spatial distributions for the introduced individuals; for each we minimize the associated cost, constrained by the requirement that the species must be restored throughout the habitat. First, we consider spatially inhomogeneous, unstable stationary solutions of the model’s equations as plausible candidates for small restoration cost. Second, we use numerical simulations to find the smallest rectangular cluster, enclosing a spatially homogeneous population density, that minimizes the cost of assured restoration. Finally, by employing simulated annealing, we minimize restoration cost among all possible initial spatial distributions of females and males. For biased sex ratios, or for a significant between-sex difference in mortality, we find that sex-specific spatial distributions minimize the cost. But as long as the sex ratio maximizes the local equilibrium density for given mortality rates, a common homogeneous distribution for both sexes that spans a critical distance yields a similarly low cost.


Introduction
Ecological restoration aims to replenish an ecosystem's biodiversity, often responding to human-induced losses of indigenous species [1,2].When ecosystem managers reintroduce a species to its former habitat, the restoration effort's success is ordinarily defined by combined ecological and economic criteria [3].Similarly, optimizing biological-control programs may integrate impact on the target species with costs of deploying the control agent [4].
Consider an example where restoration failed.Historical records indicate that Canada lynx (Lynx canadensis) were found in New York State (NYS), but were seen only rarely during most of the 20th century [7].Between 1989 and1992, no fewer than 80 lynx were captured in Canada and released in the Adirondack Mountains of NYS.Each animal carried a radio-collar, so that survival and dispersal could be monitored.The lynx rapidly dispersed; mortality during dispersal was high.Lynx population density grew too low for successful reproduction, and the species is now considered extirpated in NYS [7].
Generalizing the example, we envision restoration of a single species whose population dynamics depends on the density of each sex.Before evaluating costs, we must identify those spatial distributions of the initial population that assure successful restoration.Suppose that we initiate restoration with a single spatial cluster, within which individuals are distributed at a uniform density.Then we must find the "critical cluster" size, the minimal area the species must occupy to sustain positive population growth.Analysis of the critical-cluster criterion has advanced understanding of spatial systems in both physics [8][9][10] and ecology [11][12][13][14].However, if ecosystem managers can vary initial densities according to location, non-uniform spatial distributions might reduce restoration cost.Given multiple initial population distributions assuring sustained population increase, the most preferred option should minimize cost.Our study investigates how the minimum cost of successful restoration depends directly on spatial pattern, and how the optimal pattern depends on sex ratio, and on sex-specific mortality rates.
In this context, we model a species' restoration as a spatially detailed, two-sex reaction-diffusion system.We optimize the initial densities and spatial distributions of the sexes to minimize the cost of restoring the species to its positive, stable (homogeneous) equilibrium density throughout a habitat.The prototype of such models is the Fisher-Kolmogorov equation [15,16] which describes the dynamics of single-sex populations with logistic growth and diffusive dispersal.Our model extends the basic reaction-diffusion framework to include sex-structured dynamics [17][18][19], where an Allee effect [20,21] generates an unstable fixed point between extinction and the habitat's carrying capacity.
In many natural and managed populations, per-capitum growth is reduced as density becomes small; this is termed an Allee effect [5,22,23].Different behavioral, ecological and genetic mechanisms can induce an Allee effect [24].Low population density may diminish individual reproduction by reducing mate encounters, making prey capture more difficult, or by leaving individuals more susceptible to their own predators [25,26].Allee effects will be amplified if dispersal into unoccupied habitat reduces local population density; a sufficiently high dispersal rate can generate negative population growth, thwarting restoration [20].
Generically, the cost of restoration can be can defined as: where Ω represents the extent of the habitat, and u(r, t) is the population density at location r.Without loss of generality, we can consider a constant per-capitum cost (c * = 1).Mathematically, the restoration cost, which we seek to minimize, is a functional of the initial population's spatial distribution.The constraint requiring population persistence cannot be expressed analytically, since that would require solving the model's partial differential equations exactly.Therefore, we develop a population dynamics with simple processes and easily interpreted parameters, and obtain the minimum cost with analytical and numerical techniques.We organize the rest of the paper as follows.First, we introduce our sex-structured population dynamics, and outline the analytic and numerical methods we employ in this paper.In the Results section we conduct a systematic study of the restoration cost in multiple stages.Starting with a simplified, single-sex version of our model, we derive an unstable, aperiodic stationary solution for the PDE, resulting in a single-sex cost.We then continue with the sex-structured model and analyze the way restoration cost depends on model parameters, given a simple, constrained (i.e.homogeneous) initial distribution of individuals.Thereafter, we relax the constraints in two steps, and study how the cost can be reduced as a result.Finally, in the Discussion we compare the minimum costs found in each stage, and conclude which approach yields the most economical restoration.

General assumptions
Our model's key parameters include the sex ratio at birth and sex-specific mortality rates.We assume that females and males disperse independently by homogeneous diffusion, and that males encounter females as a mass-action process, equivalent to random mating [27].The fraction of matings leading to successful reproduction is proportional to the unoccupied fraction of the environment, (1 − m − f ) [17][18][19].That is, the population grows in a self-regulated manner.Hence, we have: where f (x, t) and m(x, t) denote the local densities of females and males, respectively.Diffusion rates are described by coefficient D f for females and D m for males.Three parameters characterize the local dynamics: θ, µ m , and µ f , which denote the fraction of individuals born female, and density-independent mortality rates for males and females, respectively.Note that while Eqs.(3) do incorporate spatial effects, they retain a "mean-field" character (the statistical physics terminology) in that correlation functions of the underlying stochastic individual-based model are factorized into products of densities [28,29].
(Note, however, that the rigorous derivation of such stochastic partial differential equations can be rather challenging [37][38][39].)Our model [Eq.( 3)], with spatial detail and diffusion removed, is identical to that studied by Tainaka et al. [17].In our earlier work [19], we briefly analyzed the above two-sex dynamics with diffusive dispersal, and found the critical radius (hence, critical cluster size) for an initial population's successful invasion of a two-dimensional habitat.While our reaction-diffusion model includes numerous simplifications for detailed application to particular species, nevertheless it exhibits the essential ecological characteristics of more complex two-sex models.Hence, implications of our results for restoration will likely hold across a wide range of specific models.
We can transform Eq. ( 3) to interpret it as a single-sex model by making the equations symmetric.To do so, we must let θ = 0.5, restrict µ f = µ m =: µ, D f = D m =: D and use the same initial-density distributions for both males and females.In this way, the two densities behave identically over time: f (x, t) = m(x, t) = u(x, t), described by the following equation: This transformation bridges the single-sex and two-sex models; using the constants in Eq. ( 4), we can directly compare results between models without rescaling parameters.Note that the (cubic) local dynamics also retains an Allee effect.

Analytic and numerical methods
Our first approach to cost minimization selects suitable unstable stationary solutions of the PDE model as initial population distributions, since we can derive them analytically for the single-sex model.These are special "critical" solutions, which transform to stable equilibrium solutions (either persistence or extinction) depending on a small perturbation.The cost associated with each stationary solution is found by numerical integration of the density profile (i.e., the area under the curve).We obtain the stationary solutions for the single-sex model by setting the left-hand side of Eq. (4) to zero, and deriving a relationship between the stationary solution's local value and its derivative.
For the sex-structured model, we cannot derive stationary solutions analytically.Instead, we base our study on the model's critical limits of cluster size and density, which are directly related to the minimum cost of successful restoration.For a homogeneous initial spatial distribution with a specific population density, there exists a critical cluster size, the smallest spatial extent such that the given density achieves sustained positive growth.Symmetrically, for a given cluster size, there is a critical initial density, the lowest density assuring sustained population growth.In both cases, the critical limit also corresponds to the minimum cost, since the cost is proportional to both cluster size and density.
The exact, parameter-dependent values of these critical limits cannot be derived analytically.Instead, we use binary search across a range of possible values, accomplished by testing each value for successful restoration vs extinction.
We discriminate restoration from extinction by numerically integrating the model until it has converged to a global equilibrium.We use second-order finite difference discretization for spatial derivatives and explicit Euler method for integration over time [40], with a sufficiently small time-step.Integration stopped when all time-derivatives at all spatial coordinates were less than 10 −6 .Finally, the cost is obtained by multiplying the initial density and the initial cluster size, and summed for males and females.
Next, we extend the uniform-density approach by allowing the initial cluster size to differ between males and females.In this case, we first calculate the critical cluster size for males at a fixed cluster size for females using binary search, resulting in a cost with respect to the given female cluster size.Then, we use gradient descent on this cost function to minimize it with respect to female cluster size.Note, that during the gradient descent we always change the female cluster size by one unit of spatial grid resolution, and we keep moving toward the negative gradient even if the local derivative is zero, because a small slope discretized with any grid can result in zero local gradients before reaching the actual minimum point.Note also, that we strongly rely on the convexity of the cost function with respect to male and female cluster sizes.
Finally, we relax all constraints on the shape of the initial population distribution; we optimize the spatially discretized shape for lowest cost under the constraint of successful restoration.Discretization is essential because it allows us to express the cost as an n-dimensional function, instead of a functional, where n is the size of the spatial grid.We use the same grid for shape discretization and numerical PDE integration, for practical reasons.Given the discretization, we use simulated annealing [43] for optimizing the shape function.This is essentially a Monte Carlo simulation, where random changes of the initial spatial distribution (the shape function) are accepted or rejected according to a specified acceptance probability function, such that the visited cost states have a Boltzman-distribution characterized by a temperature-like parameter.As this parameter is lowered, the expected value of the cost is also lowered, eventually leading to the globally optimal, minimum cost state.For the specific steps of simulated annealing, see Supporting Information (Section S5).In order to determine whether the constraint of successful restoration is satisfied we must numerically integrate the model at every Monte Carlo step.To reduce computational time, we accelerated the PDE integration with GPGPU computation using CUDA [41,42], which locates the global equilibrium of the PDE in a fraction of a second, giving a total time for the simulated annealing in the order of a few hours.

Results
The prerequisite for analyzing the cost of restoration is to ensure the local stability of successful restoration, i.e., a positive stable fixed point of local dynamics.The necessary stability condition for the two-sex model [Eq.( 3)] is given in our previous work [19]: Similarily, we can find the necessary stability condition for the single-sex model [Eq.( 4)]: These conditions provide us guidelines for selecting proper model parameters when evaluating costs in the following sections.Justification for Eqs. ( 5) and ( 6), and formulas for the fixed points are presented in the Supporting Information (Sections S1 and S2).

Unstable stationary solutions
An unstable stationary solution of the PDE seems a good candidate for the initial spatial distribution.
It is a critical solution, in the sense that given a small perturbation, the unstable stationary solution transforms to a stable, spatially homogeneous solution.Positive perturbations result in positive homogeneous population densities (successful restoration), and negative perturbations result in zero densities (extinction).Partial differential equations, such as Eqs.( 3) and ( 4) have infinitely many unstable stationary solutions.We cannot find these solutions directly, but we can derive analytical formulas for the relationship between the density and its spatial derivative.For the single-sex model [Eq.( 4)] we have the following definitions: Using v, we change variables to write a first-order differential equation: By separating variables, we obtain the following analytical solution: where E is a free parameter.The phase diagram for this equation is depicted in Fig. 1(a), and its contents are summarized as follows.
The fixed points P i correspond to homogeneous stationary solutions.Naturally, these are also fixed points of the original equations [Eq.( 4)], but here they are only special cases of stationary solutions that do not vary spatially.Hence v(u) = 0. P 1 and P 3 are saddle points (stable equilibrium nodes of the local dynamics) corresponding to extinction (u = 0) and persistence (u > 0), respectively.P 2 is a center (unstable fixed point of the local dynamics) corresponding to the unstable equilibrium due to the local dynamics' strong Allee effect.
Curves in Fig. 1(a) correspond to inhomogeneous stationary solutions that can be classified by the value of the free parameter E in Eq. (10).Separatrices S 1 , S 2 , and S 3 correspond to the following values (with the same subscripts, see Supporting information Section S3 for details): where u * 1 , u * 2 , and u * 3 are the equilibrium densities of local dynamics, i.e., the values of u corresponding to P 1 , P 2 , and P 3 , respectively.The closed elliptical curves around P 2 represent periodic stationary solutions, and they are the only ones of interest, because all other curves extend to infinitely large negative or positive densities, neither having biological meaning.
Spatially periodic stationary solutions may offer candidate initial population distributions.In principle, if minimum densities within each period were close to zero, then we could select a segment of the solution, one period in length between two density minima, and apply it as an initial spatial distribution.However, in our case, as the minimum value of u goes to zero, the period of the solutions goes to infinity, and the curves converge to the S 2 separatrix, which corresponds to an aperiodic stationary solution.The exact aperiodic shape of u(x) can be found by numerical integration along S 2 , depicted in Fig. 1(b).Since u(x) converges to zero rapidly, we can use it as an initial population distribution by taking its central segment above an arbitrary small (biologically meaningful) density threshold.Because of the fast convergence to zero, the length of the segment will be finite.In our study, we use 10 −5 for the density threshold.Then, for every D and µ combination, we have an exact shape, and an exact cost value defined as twice the area under the curve (counting both males and females), denoted by C stat .We will compare these costs with those found by the other two methods described in the following subsections.
As an interesting observation, we note that contrary to intuition, the period length of the stationary solution does not converge to zero as the solution curves approach P 2 , see Supporting Information (Section S4).
For the sex-structured model [Eq.3] we cannot derive unstable stationary solutions analytically.Instead, we study the scaling of restoration cost through a numerical analysis of the critical cluster size.

Critical cluster size and minimum cost
Criticality of an initial cluster's size occurs when density reduction due to dispersal exactly balances the net effect of local natality and mortality; Supporting Video S1 clearly shows this form of criticality in one dimension.In two dimensions the expanding population front's speed is reduced in proportion to the curvature of the cluster; therefore, it affects the critical cluster size.To avoid confusing effects of curvature with other parameters' impact on the critical cluster size and, hence, the minimum cost, we restrict our study to one-dimensional space.
We begin our analysis of the critical cluster size by assuming a uniform spatial distribution.This is the most obvious choice, for its mathematical simplicity, and its plausible application (e.g., an animal population surrounded by a fence before release can be modeled with a uniform spatial distribution).Therefore, we have four parameters describing the distribution: f 0 , m 0 , l f , l m , which represent female density, male density, length of space occupied by females, and length of space occupied by males, respectively, at the initiation of restoration.We refer to this initial setup as "rectangular", for its shape on a density vs location plot, and we shall refer to l f and l m as the cluster sizes of the initial population.The cost of the initial state is defined simply as: and it is minimized by using the critical cluster sizes l * m and l * m , for males and females, respectively.Note, that l * m and l * m are themselves dependent on initial densities f 0 and m 0 , as well as model parameters D f , D m , µ f , µ m , and θ.Therefore, we systematically study dependence of critical clusters on all parameters.
As an initial step, we analyze the dependence of critical cluster size on the diffusion coefficients.We anticipate that the critical cluster size, i.e., critical length (in one-dimensional space) is proportional to the square root of the diffusion coefficient.Similar scaling has been observed in two-dimensions when a population with a strong Allee effect disperses by diffusion [19,20].Our aim here is to show the same behavior in one-dimensional space, and to ask whether it holds when male and female diffusion coefficients and cluster sizes differ.
Figure 2(a,b) shows that as long as we employ the same cluster size for males and females, the critical cluster size has the expected scaling behavior l * ∼ √ D [20] with respect to both male and female diffusion coefficients, even if one sex has a fixed diffusion coefficient.Note that it is sufficient to study the dependence on the diffusion of one sex (here, males) while the other is fixed (here, females), because of the symmetric construction of the model.Further, Fig. 2(c,d) indicates that fixing the cluster size of one sex while measuring the critical cluster size of the other sex with respect to diffusion coefficients yields non-trivial scaling.However, we observe that a higher dispersal rate results in a larger critical cluster size, and, hence, a larger restoration cost.
Continuing our analysis of critical cluster sizes, we now assume that the critical cluster sizes are equal for both sexes (l * f = l * m =: l * ); we relax this constraint later.Even with the equal cluster-size constraint, initial densities for males and females within that cluster may, in principle, differ.In practice, a density difference could be implemented for most diecious species.To find the best choice of initial density values (with respect to minimizing cost), we aim to relate them to model parameters, taken as given for the focal population.
Naturally, the initial densities must exceed the Allee threshold (the unstable fixed-point densities); otherwise, the population can never achieve positive growth.Figure 3(a) shows that as the initial density is lowered, the critical cluster size increases, and goes to infinity as we approach the Allee threshold.
Scaling of the cost, however, is non-trivial.We can always find the minimum cost at a density somewhere in the vicinity of the positive stable fixed point of the system, but always slightly below it.We understand this by considering the dynamics just after initial introduction.If the population starts from its stationary density (the stable, positive fixed point), then the local densities can only decrease, due to diffusive dispersal.However, if the initial density is lowered slightly, then the population has a chance to grow locally (in particular, at the center of the cluster) before the effects of diffusion reach it, while the eventual spread through the habitat remains the same.In essence, the cost is slightly lowered by handing over some of the spreading effort to growth dynamics.However, as Fig. 3(b) and Fig. 3(c) show, this advantage in cost-reduction is very small.We can conclude that using the stationary densities as initial densities results in a sufficiently low cost.Also, it provides a good choice of initial density based on model parameters, because the stationary densities depend on local dynamics, which in turn depend on model parameters.The formulas for the stationary densities are included in the Supporting Information (Sections S1 and S2).
We expand on the cost-minimizing property of stationary densities; we use them as initial densities throughout the rest of our study.We continue with the analysis of the critical cluster size's dependence on model parameters.In particular, we focus on the value of sex ratio at birth (θ), because both population stability and equilibrium density depend on this single parameter.For population stability and persistence, there is a range of permissible sex ratios, determined by the mortality rates (see Supporting In every case, initial (spatially homogeneous) densities are fixed: f 0 = m 0 = 0.45.
Information Section S2).Within this range, we define the optimal sex ratio θ * as the value maximizing the equilibrium population density [17,19]: Note, that when the mortality rates are equal, the optimal sex ratio is 0.5, which is the parametrization, by definition, in the symmetric single-sex model.We find that the smallest cluster size assuring restoration corresponds closely to the optimal sex ratio, and that any small deviation from the optimal value causes a small increase in the critical cluster size [Fig.4(a)].However, since the equilibrium population densities (serving as initial densities) decrease at suboptimal sex ratios (by definition; see Fig. S1 in Supporting Information), the combined effect on the cost is non-trivial.As we see on Fig. 4(b), restoration cost is minimized at approximately the same sex ratios minimizing the critical cluster size, indicating that the cluster size is more sensitive to biased sex ratios than to equilibrium densities.We also find that strongly biased sex ratios approaching the boundary of the stability range cause both cluster size and restoration cost to diverge.
To complete the relationship between the sex ratio minimizing restoration cost ( θ) and the sex ratio that locally maximizes total equilibrium population density (θ * ), we compare the two quantities numerically.For the latter we have an analytical expression [Eq.(15)].But our cost-minimizing sex ratios have only limited precision, since for each value of θ we employed binary search to determine the critical cluster size, which, in turn, determines the cost.Therefore we define a computational error bound on θ as the range of θ values that give critical cluster sizes within the error range of the minimum point's cluster size found by binary search.(15)] and numerical bounds of cost-minimizing sex ratios θ, using homogeneous population distributions with stationary initial densities, D = 1.0,(a) µ m = 0.01, (b) µ m = 0.02, (c) µ m = 0.03.
Figure 5 offers comparison of the density-maximizing and cost-minimizing sex ratios.We conclude that the sex ratio maximizing equilibrium density is identical to the sex ratio minimizing restoration cost, up to computational error.
To this point, our results reflect the assumption that individuals of each sex are introduced across the same extent of habitat.We now relax this constraint; that is, we permit l * f = l * m , and ask whether the cost of restoration can be reduced by introducing individuals into sex-specific lengths of habitat.We denote the ratio of costs obtained by unequal and equal cluster sizes as: Figuse 6 presents our results.We find that optimality of the sex ratio plays a central role in reducing the cost of restoration.In particular, when the sex ratio equals the optimal (density-maximizing) value determined by mortality rates, the minimum cost is achieved with equally sized clusters.No further cost reduction can be achieved by allowing different cluster sizes for males and females.However, the relative advantage of sex-specific cluster size increases as the sex ratio deviates from the optimal value.If we consider that the absolute cost value diverges for strongly biased sex ratios [see Fig. 4(b)] we conclude that in such cases the savings achieved by adjusting the initial cluster sizes of each sex could be substantial.

Simulated annealing
We now relax all constraints on the initial cluster's spatial distribution, and use simulated annealing to optimize sex-specific distribution shapes with respect to cost.At this stage we assume only that the costminimizing distributions have a finite support, and we carry out the minimization accordingly.However, the support of the function is allowed to grow or shrink by random shape changes during simulated annealing; see Supporting Information (Section S5) for details.By analyzing a series of minimum cost distributions obtained with simulated annealing, we observe the following.First, the distributions indeed have a finite support.Although we initialize them as such, the width of the optimal population distribution tends to become smaller, rather than larger, during simulated annealing.This effect during the minimization procedure is shown by Supporting Video S2; typical final, optimized shapes are shown in Fig. 7.It is also remarkable that the edges of the distributions go to zero very sharply; this property develops without any influence inherent to the procedure.
Generally, the final result is an "arch"-shaped distribution, with similar dimensions for females and males.Note that as the sex ratio diverges from its optimal value, we observe a change in the sizes of the two initial population distributions and in the height of the peaks.These changes in spatial distributions occur roughly in proportion to the system's positive stationary densities [Fig.7(b,c)].Interestingly, the height of the peaks always falls between the stationary densities and the Allee threshold.Note that this shape provides the maximal rate of population growth possible during the first moments of the simulation, hence it combats the diffusion-amplified Allee effect most efficiently.Costs corresponding to the optimized distributions are denoted C SA when we compare results with other methods.

Discussion
We examined three approaches to minimizing the cost of a species' restoration; the approaches differ in both ecological premises and mathematical methods.We considered the aperiodic, spatially inhomogeneous solution to the single-sex dynamics [Fig.Summary comparisons of the minimum restoration cost achieved by the different methods for the symmetric single sex model appear in Fig. 8.The aperiodic stationary solution to the dynamics gives significantly larger cost than the other approaches.Considering the shape of these distributions, they would likely prove difficult to implement in application.Restricting attention to the other two approaches, it is remarkable that simple homogeneous distributions and the results of simulated annealing yield essentially identical costs.
Figure 9 compares minimum costs for each critical-cluster analysis (i.e., a single cluster size and sex-specific cluster sizes) and costs incurred under simulated annealing.The critical cluster methods assume uniform initial population density within cluster bounds; simulated annealing lets initial densities depend on spatial location.The minimum cost varies little among methods as long as the sex ratio at birth does not deviate too much from the optimal value (here, θ * = 0.5).As sex-ratio bias increases, optimal sex-specific initial cluster sizes can lower the minimum cost of restoration.Simulated annealing reduces restoration cost even further, but this advantage becomes significant only at strongly biased (and biologically rare) sex ratios, and implementing such spatial distributions in application could prove difficult, negating any cost advantage.The same qualitative conclusions hold when we fix the sex ratio and increase the difference between the sexes' respective mortality rates, because the mortality bias can also increase the difference between the optimal and any fixed sex ratio.
Our model assumes deterministic dynamics, which does not account for extinction due to demographic stochasticity in populations near an extinction threshold [44,45].This effect can be exaggerated when a population's spatial dispersion leaves dynamically independent clusters near critical size [12,46].
We assume diffusive dispersal.Many plants, and some animals, are "diffusion limited," i.e., the probability of long-distance dispersal is much lower than diffusion assumes [47][48][49].Diffusion limitation implies that the criterion for expected positive growth demands greater propagation, relative to mortality, as dispersal limitation increases [13,50].We also assume that no explicit interspecific interactions affect the population during restoration.Species occupying the community to be restored may facilitate restoration; for example, trees may attract birds that disperse seeds of other tree species [51].Alternatively, resident species may resist the introduced species biotically [52,53].Interspecific interactions will often affect the likelihood of restoration success, as well as the cost.Consequences of these interactions can sometimes be expressed abstractly through the introduced species' positive equilibrium density; in other cases, successful restoration may demand quantification of these interactions.
We assume an Allee effect arises from interaction of self-regulation with a birth rate that depends on the density of each sex.In the context of restoration, a two-sex dynamics may be essential to predicting spatial-expansion rate if dispersal differs between sexes [6].We model mating encounters via massaction, which should be reasonable for animals maintaining individual home ranges, or for dioecious plants with random mating.Alternative "marriage functions" [54] apply to certain species, particularly for polygynous or polyandrous mating systems.
We modeled a single species' restoration only.Habitat restoration may attempt to manage particular multi-species interactions, or may seek to promote growth of many threatened species [1].Our cost function ignores feedback of a species' restoration on other biotic processes, or on economic stake-holders who incur post-restoration costs [55].
Our results suggest some considerations for species restoration.First, if a species disperses rapidly, individuals should be introduced concurrently, rather than serially.The initial population will increase only if the density exceeds any Allee threshold, and continues to do so as individuals disperse.Intuitively, the number/density of individuals introduced should increase with their dispersal rate.
Second, restoration cost declines little by introducing a species at a density below the (estimated) carrying capacity, unless the species disperses very slowly.Third, a uniform spatial distribution within the initial populations's cluster adds little or no proportional cost over the ogive profile assumed in our simulated annealing method, as long as the sex ratio is close to optimal.Spatial uniformity will likely prove more practical for most animals.Given a uniform density close to the positive, stable equilibrium, restoration should focus on an initial population whose expanse exceeds the critical-cluster size, which (again) increases with dispersal rate.
Finally, adjusting frequencies of the sexes in an initial population may decrease the cost of successful restoration.Of course, if one sex always limits population growth, an excess of that sex promotes restoration.If population growth depends on the density of each sex, introducing the sexes at different densities, or with different cluster sizes, may prove advantageous.Sex ratio at birth may be unbiased, but mortality rates may differ between sexes, particularly during dispersal.Adjusting the sex ratio at introduction to match frequencies at positive equilibrium densities should promote successful restoration, and reduce its cost.

S4 Analyzing stationary solutions of the single-sex model
Stationary solutions outsite S 2 separatrix (blue, green, and black curves on Fig. S2) have no physical or biological meaning, because they are unbounded, and extend to negative density ranges.Therefore we restrict our analysis to the periodic stationary solutions (red curves), and show that there is a minimum required spatial extent for their existence.First, we find the stability matrix and its eigenvalues.We obtain the following eigenvalues: We are interested in the imaginary eigenvalues corresponding to the center (P 2 ): As we approach P 2 , the orbits around it converge to a circle, with the angular frequency given by the imaginary eigenvalue: The periodic solution has the same angular frequency, therefore the wavelength in the limit of P 2 is: This is the minimum habitat size allowing periodic stationary solutions.Note, that this case corresponds to E = E 1 .

S5 Simulated annealing
To move beyond the homogeneous distribution of the initial population as we seek to minimize cost, we need a method that finds the minimum cost in the domain of infinitely many possible distribution functions.The key to solving this problem is to define the initial distribution in discretized form.We can determine the eventual persistence or extinction of an initial population only by running the simulation to global equilibrium, and this simulation requires a spatial discretization for the integration.We use the same grid to discretize the initial distribution function.Further, we can assume that whatever initial distribution gives the minimum cost, it will have a finite support, therefore we can restrict the spatial extent of the distribution to a region in the center of the simulated space.If this restricted spatial domain includes k grid points, we essentially reduce an infinitely detailed initial distribution to a k-dimensional function.
To find the minimum cost distribution, we use a global optimization method on the k-dimensional cost function, given the constraint that the population must persist.We use simulated annealing [4], which

Figure 1 .
Figure 1.(a) Phase plot of stationary solutions in the single-sex model, described by Eq. 10.D = 1.0, µ = 0.05.The dots indicate fixed points, the thick lines indicate separatrices.Different curves correspond to different E parameters, however, values of E were not chosen uniformly, for aesthetic reasons.(b) The stationary solutions found by integrating along the S 2 separatrix, for multiple mortality rate parameters; D = 1.0; x is distance form the habitat's center.

Figure 3 .
Figure 3. Scaling of (a) critical cluster length at introduction, and (b) cost, with respect to initial population density.The total density shown on the x axis is divided equally between males and females.The point markers on the curves show the stable stationary densities of local dynamics.(c) shows the cost landscape with respect to male and female densities; the blue cross marks the stable stationary density values; µ f = 0.04, µ m = 0.03, θ = 0.5.For all figures, D = 1.0.

Figure 4 .Figure 5 .
Figure 4. Scaling of (a) critical cluster size and (b) cost, with respect to sex ratio, at different mortality rate combinations; D = 1.0.Rectangular initial populations were used with stationary population densities.

Figure 6 .
Figure 6.Minimum cost found by allowing different initial cluster sizes for males and females, relative to the case where cluster lengths are equal.Common parameter: D = 1.0.Individual parameters: (a) µ f = 0.01, (b) µ f = 0.03.