Modelling three-dimensional space to design prey refuges using video game software

. Refuges can be ecologically important, allowing access only to some species or individuals and providing prey protection from predators. Creation of refuges can be used to protect threatened species from introduced predators, which can have large negative impacts that are dif ﬁ cult to attenuate via other means. To design refuges for conservation purposes, refuge accessibility to different species must be understood. Traditional techniques are not adequate to measure or describe complex three-dimensional spaces which are often important refuges. We designed a novel predictive method for modeling three-dimensional refuge space using video game software that simulates real-world physics (Unity, PhysX). We use the study system of endemic New Zealand skinks ( Oligosoma spp.), their introduced predators, house mice ( Mus musculus ), and the habitat of interstitial spaces within rock piles to demonstrate how this modeling technique can be used to inform design of habitat enhancement for conservation. We used video game software to model realistic rock piles and measure their interstitial spaces, and found that the spaces we predicted matched those we measured in real rock piles using computed tomography (CT) scanning. We used information about the sizes of gaps accessible to skinks and mice and the results of our modeling to determine the optimal size of rocks to create refuges which would protect skinks from mice. We determined the ideal rock size to be those with graded diameters of 20 – 40 mm. The approach we developed could be used to describe interstitial spaces in habitats as they naturally occur, or it could be applied to design habitats to bene ﬁ t particular species.


INTRODUCTION
Introduced predators can have detrimental effects including causing the extinctions of many species, and pose a challenge for wildlife managers seeking to protect threatened endemics (Clavero and Garc ıa-Berthou 2005). Introduced predators are extremely difficult to control or eradicate over large areas, meaning that managers must look to other strategies to reduce their negative impacts (Hulme 2006). The presence of refuges that can be used as shelter can reduce the detrimental effects of predators, including those of introduced predators, on prey (Stuart-Smith et al. 2007). Optimal foraging theory predicts that predation depends on the costbenefit ratio involved in acquiring and consuming prey (Pyke et al. 1977); the addition of refuges may increase the cost of predation for predators and therefore reduce its occurrence (Whitlow et al. 2003). Creation of refuges via habitat manipulation thus represents an opportunity for conservationists seeking to protect threatened species (Sinclair et al. 1998), and has been used for the conservation of many taxa including reptiles (Croak et al. 2013), mammals (Rouco et al. 2011), birds (Belthoff and Smith 2003), invertebrates (Green 2005), and aquatic species (London Convention and Protocol/UNEP 2009). Previous approaches for measuring and characterizing refuges, however, are insufficient, particularly for species requiring complex threedimensional structures.
Approaches to characterizing refuge space have previously fallen into three categories. The first is direct physical measurement of refuges in situ, for instance, using probes to measure the length of lizard burrows (Milne and Bull 2000); expanding foam casts to measure spaces under rocks (Croak et al. 2008); or grids of gauges to measure reef topography (Alexander 2013). This approach allows researchers to characterize natural refuges, and can be useful for correlative investigations of refuge use. However, physically measuring three-dimensional space can be difficult and requires either relatively simple refuges (e.g., a cylindrical burrow) or the use of simplified measures (e.g., gauges measure reef topology from above, but cannot reach any spaces below overhangs). A second approach is to create artificial refuges, the shape, and size of which can be experimentally varied (Warfe et al. 2008, Hesterberg et al. 2017). This approach gives researchers more control over variables, allowing for more rigorous experiments. To allow for clear comparisons of size and/or shape, however, these refuges tend to be simplified. Studies using this approach thus tend to report the effects of artificial refuges on animals rather than describing the refuges themselves. A third approach is to characterize habitat complexity by using a ratio of refuge size to predator size to predict maneuverability of the predator and therefore the prey survivorship (Bartholomew et al. 2000). The use of the average size of spaces is a major limitation of this approach, as variably arranged spaces can result in the same score despite having vastly different effects on predator maneuverability (Bell et al. 2003).
Recently advances have been made in understanding three-dimensional refuges in aquatic systems, for instance the use of computed tomography (CT) scanning or photogrammetry to generate three-dimensional models of mangrove roots (Kamal et al. 2014), and kelp holdfasts (Orland et al. 2016). Ware et al. (2019) used spherical space analysis to estimate accessible and inaccessible refuge space within macroalgae based on a three-dimensional computer model of macroalgae. Of these approaches, computerbased modeling is preferable as it is cheaper than CT scanning, and has the potential to be used more widely in conservation management.
New Zealand is home to a large diversity of endemic lizards from two families, geckos (Diplodactylidae) and skinks (Scincidae), which are highly threatened (Hitchmough et al. 2016a). New Zealand lizards are particularly vulnerable to predation by introduced mammals due to their evolutionary naivety (Hitchmough et al. 2016b). House mice (Mus musculus) are introduced predators of skinks in New Zealand (Newman 1994, Wedding 2007, Norbury et al. 2014. Mouse population explosions have been responsible for, or correlated with, declines in populations of several skink species (Newman 1994, Ussher 2006as reported in Wedding 2007. Mice are harder to control than larger invasive mammals and are therefore seldom targeted in control operations (Innes et al. 2012, Hitchmough et al. 2016b. Conserving endemic lizards in New Zealand will thus require alternative management strategies aimed at mitigating the effects of mice. Many endemic New Zealand skink species live in complex three-dimensional rocky habitats, where they bask on exposed surfaces and shelter in interstitial spaces (Towns 1991, 1996, Towns and Ferreira 2001, Hoare et al. 2007). Skinks' natural predator avoidance behavior of sheltering inside rock piles can allow them to escape predators such as birds and larger mammals; however, it may do little to protect them from mice, which are small enough to enter many of the same interstitial spaces as skinks and have been directly observed to attack skinks basking near rock refuges (Norbury et al. 2014).
v www.esajournals.org Management of skinks in New Zealand sometimes involves building new rock pile habitat (Anderson et al. 2012), which is currently based on best guesses by managers, with no scientific testing of rock pile composition. The characteristics of refuges in the new rock pile habitat may affect skink/mouse interactions and skink survival. Creating rock piles that are accessible to skinks but not to mice would provide natural protection for skinks by providing refuges that allow skinks to avoid mice, and by hindering the movement and occupancy of mice within the piles, reducing the likelihood of creating additional mouse habitat (particularly for nesting). This would increase the cost-benefit ratio of mouse predation on skinks, reducing opportunistic predation by mice, and increasing skink survival, recruitment, and population viability. However, the interstitial spaces within the rock piles cannot be measured using traditional means.
We present a novel computer-based approach for generating realistic rock pile models and measuring ecologically relevant aspects of threedimensional refuge space. Our method uses physics engine software designed for video games (NVIDIA Corporation 2011, Unity Technologies 2016) to predict the sizes of interstitial spaces within rock piles, but could be generalized to characterize other kinds of refuge space. We test the accuracy of our predictions using CT scanning of rock piles. The results can be used to inform refuge design, which we demonstrate by using this technique to design refuges to protect endemic New Zealand skinks (Oligosoma spp.) from predation by introduced mice. As mice are slightly larger than many New Zealand skink species, we aim to design rock piles which optimize interstitial spaces that are accessible to skinks but are too small for mice.

Sizes of gaps accessible to skinks and mice
In order to determine the accessibility of rock piles, we needed information about the sizes of gaps that skinks and mice can enter. Information for mice was available from published experiments for designing predator exclusion fencing, testing the ability of mice to escape confinement through barriers of varying sizes (Day and MacGibbon 2007). Comparable data for skinks were unavailable; we therefore performed an experiment to determine skinks' ability to escape a box through variably sized holes (Appendix S1).
For any animal, access to a space will be constrained by the animal's size. There is therefore a minimum gap size below which the animal will be unable to access a space. Where there are two size classes of animal, gaps can be either too small for either animal to access, large enough for the smaller animal but too small for the larger animal, or large enough for either animal to access. We named these categories small, optimal, and large, respectively (optimal refers to access for the smaller [prey] animal; in this example, skinks). Classifying gaps into these three categories enables an understanding of their accessibility for each animal. We used gap accessibility data for skinks and mice to define the three size categories as follows: small, any dimension ≤3 mm; optimal, all dimensions >3 mm, at least one dimension <10 mm; and large, all dimensions ≥10 mm (Appendix S1).

Virtual modeling of three-dimensional refuge space
Summary of methods.-We used the free video games design program, Unity (Version 5.5.5f1; Unity Technologies 2016), and its inbuilt physics engine PhysX (Version 3.3 SDK, NVIDIA Corporation 2011) to model the behavior of virtual rocks dropped from a height and allowed to settle into a pile. Physics engines are simulation software which model realistic movement of objects by simulating physical traits such as mass, gravity, friction, inertia, momentum, and collisions between dynamic objects. We then used a technique called ray casting, also within Unity, to measure the interstitial spaces in the simulated rock pile. We simulated rock piles of varying compositions (sizes, shapes, and combinations of rocks). We used R (version 3.5.3; R Core Team 2016) to translate the raw data obtained from Unity into three-dimensional measurements for spaces inside the simulated rock piles. We categorized the measured spaces into small, optimum, or large size classes, and calculated the proportion of spaces in each pile that fell into each size category, to estimate the v www.esajournals.org accessibility of rock piles to skinks and mice. To check that the results from our modeling were accurate for real rock piles, we CT-scanned rock piles and digitally measured the three-dimensional images using ImageJ (version 1.50i; Schneider et al. 2012), and compared the results with the predictions from our modeling.
Rock models.-To simulate physical objects, PhysX and Unity use rigidbodies which are virtual three-dimensional objects that can receive forces and torque to act in a realistic way. PhysX models rigidbodies as either predetermined simple shapes or more complex meshes. We used the free software Blender (version 2.77; Blender Foundation 2016) to make meshes that were modeled after real rocks collected from a quarry ( Fig. 1). We made 20 rock meshes, and we stretched/squashed them in three dimensions to vary their sizes and shapes. Deforming the rocks' shapes allowed us to have a larger variety of individual rock models, and changing their sizes allowed us to test different grades of rock. Each mesh became a rigidbody in Unity, resulting in rock models that behaved naturally according to physical forces. PhysX does not allow individual rigidbodies to be concave; because of this, we modeled the rocks as convex shapes only.
The size of each rock model is represented as a range: This is the equivalent of the grading used to measure rocks in construction and landscaping (Table 1). The two numbers indicate the sizes of two screens that can be used to sort the rocks, for example, a rock graded 20-40 mm will go through a 40-mm screen but not a 20-mm screen.
Dropping the rocks.-We wrote a C# script for Unity (Appendix S2) which simulated tipping a number of rocks from a height and allowing them to settle in a pile before they are measured using ray casts ( Figure 1). The script allowed us to control aspects such as the rocks' sizes and shapes, the number of rocks, the combination of different types of rocks, and the height from which they are dropped.
We tested 14 different rock compositions ( Table 1). The number of rocks per pile varied from 3500 to 6000 depending on their sizes; we generated more rocks when the rocks were smaller in order to make piles of approximately the same size. The piles we generated were~1-2 m in diameter and 0.5 m high. We ran the program 35 times for each composition, generating and measuring 35 rock piles with the same types of rocks. This was because by 35 simulations the means of gap sizes had reached an asymptote. Following an exploratory assessment of various sizes and shapes of rocks, we performed a more exhaustive qualitative assessment of three different scenarios: (1) Rock piles comprised entirely of one grade of rock; (2) rocks of grades 11-20 and 70-100 mm in combinations of different proportions; and (3) rocks of grades 10-14 and 20-40 mm in combinations of different proportions (Table 1). We chose these grades based on outcomes from the pilot study indicating the approximate appropriate size range, and because 10-14 and 20-40 mm grades are commercially available. We also tested three combinations in order to compare them with CT-scanned exemplar rock piles (Table 1).
Measuring gaps: ray casting and analysis.-We wrote a C# script for Unity which uses ray casts to measure gaps within the pile in three dimensions (Appendix S2). Ray casts are virtual rays that have a defined point of origin, direction, and length. Each ray cast records the position of every point along its length at which it intersects the edge of an object. The script generates three groups of ray casts which penetrate the virtual rock pile in three orthogonal directions: an X group, a Y group, and a Z group. Each group's rays' points of origin lie on a two-dimensional grid oriented perpendicular to the direction the rays are traveling in. This arrangement results in a three-dimensional grid of points within the rock pile at which rays from all three directions intersect ( Figure 1). Some points fall within interstitial spaces and some inside rocks; we used R to discriminate between these so we only analyzed points that fell within interstitial spaces. This ray casting method means that individual interstitial spaces may be sampled more than once, and larger interstitial spaces are more likely to be sampled and to be sampled more frequently. The data output from this script is a group of files with coordinates of every point where a ray cast intersected the edge of a rock (hit locations). We used R (R Core Team 2016) to extract the X, Y, and Z sizes of all sampled gaps from these data (Appendix S3). Briefly, the R script does the following. First, it arranges the hit locations using the coordinate information so that all hit locations for each single ray are arranged in order. Second, it pairs hit locations v www.esajournals.org that bound gaps along each ray to get locations and one-dimensional sizes for each gap along each ray. Third, it compares the location data from the X, Y, and Z rays to match the three onedimensional size measurements of each sampled gap. The results are X, Y, and Z size measurements for the space around each point from the sampling grid that fell within an interstitial space.
For each pile generated, we used ray casts to sample 4000 points spaced in a grid. The grid was arranged so there were 50-mm spaces between each point. The grid started 10 mm off the ground and was 500 mm high and 1 m in width and depth, and was positioned approximately around the center of the rock pile. We chose this spacing to ensure that it was spread throughout the entire pile so that the sampling was not biased (spaces in the center may be different from those at the edges or top), and to create a moderate amount of data to work with in R (larger datasets increase the analysis time). We classified each measured gap as either small, optimum, or large using R (Appendix S3). We calculated the mean proportions of small, optimal, and large gaps for each rock pile composition and calculated 95% confidence intervals using the percentile method.

CT scanning
To compare the predictions from our modeling experiment with physical measurements, we used CT scanning to measure the interstitial spaces of real rock piles. We tipped rocks into 550 9 260 9 410 mm bins which were scanned at Massey University's Veterinary Teaching Hospital using a Brilliance CT 16 slice scanner (Philips, The Netherlands). The bins were necessary to allow the rocks to be stacked into a pile but still fit through the scanner's 700 mm diameter gantry (Figure 2). We were limited to testing three different compositions due to the expense of hiring the CT scanner.
The three compositions we tested were as follows: (1) rough rocks collected from a quarry, sorted by hand using screens, with a combination of 50% (by volume) rocks graded 10-30 mm, and 50% rocks graded 50-80 mm.
(2) Rough rocks collected from a quarry, sorted by hand and using screens, and graded 50-80 mm. This was chosen to investigate whether the predictions from the computer modeling would be accurate when rocks (and gaps) are larger. (3) Smooth flat oval river stones purchased from a landscaping supplier (beach flats, medium, purchased from The Goods Shed, Wellington, New Zealand). The rocks' three dimensions were~10-20 mm, 45-75 mm, and~50-100 mm (i.e., approximately in the size range of other rocks tested but with a more flat, round shape). We used these to investigate whether there were any adverse effects from modeling the rocks as entirely convex shapes (round rocks are entirely convex, unlike the quarry rocks we used as models which generally have some concavities).
We scanned each rock composition twice, using the same rocks and thoroughly mixing them between scans. The scans had a resolution Note: We modeled the rough shapes after rough rocks from a quarry and the smooth shapes after oval river rocks. † Grades refers to the size of each rock and is represented as a range: the equivalent of the grading used to measure rocks in construction or landscaping. The two numbers indicate the sizes of two screens that can be used to sort the rocks, for example, a rock of grade 20-40 mm will pass through a 40-mm screen but not a 20-mm screen.
‡ Where two grades are listed, two different sizes of rocks were mixed together; combination describes these mixes as a percentage of the rocks that were of the smaller grade. Ellipses indicate only one grade of rock was used for the pile.
§ Number of rocks is how many rocks were generated in the Unity simulation (we generated more rocks when the rocks were smaller so that the piles created were approximately equal sizes, 1-2 m wide and 0.5 m high). of 1 mm. The outputs were DICOM images, a medical imaging file format containing stacked slices of a three-dimensional image.
We used ImageJ (Schneider et al. 2012) to measure gaps in the CT images. We trimmed the images to include the entire width (410 mm) and height (260 mm) of the pile, but cut~40% of the length (leaving~300 mm), to make processing less computationally and labor-intensive. To make the measurements comparable to those obtained from the Unity simulation, we wrote a macro which displayed the DICOM images in three axes (XY, XZ, and YZ) simultaneously and calibrated a cursor to match within all three images (Fig. 2). A second macro placed a three-dimensional sampling grid of points within the images at 50-mm spacings. For each point that fell within an interstitial space, we used ImageJ's measuring tool to measure the X, Y, and Z dimensions of the space from that point. We measured gaps at~500-600 points within each pile. If several points fell within the same gap, we sampled the gap at each point. We used the same methods in R as we used on the data from Unity to categorize the gaps as small, optimal, or large.

Virtual modeling of three-dimensional refuge space
The proportions of gaps varied with the rock pile composition (Fig. 3). For piles composed of just one size of rock, there was a trend for smaller rocks to produce more small and optimal gaps and fewer large gaps, and for larger rocks to produce more large gaps and fewer small and optimal gaps (e.g., 45-80 mm rocks produced 2% small, 13% optimal, and 85% large gaps, whereas 11-20 mm rocks produced 25% small, 64% optimal, and 9% large gaps; Fig. 3A). For piles composed of mixes of two rock grades, there was a linear trend of the proportion of optimal gaps and small gaps increasing, and the proportion of large gaps decreasing, as the proportion of small rocks in the mix increased. For the piles composed of combinations of 10-14 and 20-40 mm Fig. 2. Screen capture from ImageJ (Schneider et al. 2012). Images from a CT scan representing the three dimensions of a rock pile. A cursor (yellow lines) is calibrated between the three windows and used to measurement interstitial spaces at a predetermined grid of sampling points. Insert bottom right: rocks stacked and ready to enter the CT scanner's gantry.

CT scanning
The observed proportions from the CT scans followed the same pattern as we predicted using the modeling method (Fig. 4). The observed values fell within the predicted 95% confidence interval for the 10-30 mm + 50-80 mm composition. For the pile with round rocks, the observed proportions of small gaps fell within the predicted 95% confidence interval; however, only one of the observed values for proportion of optimal gaps fell within the predicted 95% confidence interval. The observed values for the 50-80 mm composition fell outside of the 95% confidence value, by 2-3% for the proportion of Fig. 3. Proportion of measured gap space that is small, large, and optimal in rock piles simulated using PhysX and Unity. (A) Rock piles are composed of a single grade of rock only, and the rock size is represented as a range: the equivalent of the grading used to measure rocks in construction and landscaping. The two numbers indicate the sizes of two screens used to sort the rocks, that is, the rock will not pass through a screen of the smaller size but will pass through a screen of the larger size. Rock piles in (B) and (C) are composed of two grades of rock mixed together: Percentages represent how much of the mix is rocks of the smaller grade. In (C), 0% and 100% columns are piles composed of only the larger and smaller grades, respectively. Rock grades in (B) are 11-20 and 70-100 mm, and rock grades in (C) are 10-14 and 20-40 mm. See Table 1 for more details of compositions. Error bars are 95% confidence interval. N = 35 simulations per composition. Asterisk indicates hypothesized best design for skink refuges. Fig. 4. Proportions of optimal, small, and large spaces inside rock piles measured from CT scans. Rock size is represented as a range: the equivalent of the grading used to measure rocks in construction and landscaping. The two numbers indicate the sizes of two screens used to sort the rocks. For the composition which is a combination of two grades, the proportion of the two grades is 1:1 (by volume). River rocks are smooth flat oval river stones purchased from a landscaping supplier (dimensions~10-20 mm,~45-75 mm, and~50-100 mm). Columns are values predicted from a virtual modeling technique, with error bars for 95% confidence interval (N = 35 per composition). Xs are the proportions observed in CT scans (N = 2 per composition).
v www.esajournals.org small gaps and~10% for the proportion of optimal gaps.

DISCUSSION
PhysX has previously been used to design haptic feedback models for training medical professionals (Ricardez et al. 2018) and to simulate the interaction of ships and ice sheets (Lubbad and Løset 2011). The current study represents a novel use of this technology for ecological purposes.

Modeling three-dimensional space to understand refuges
Our results showed that it is possible to use three-dimensional modeling techniques to predict the sizes of interstitial spaces within rock piles, and to observe differences in the sizes of interstitial spaces between rock piles with different compositions. CT scanning, limited to two for each rock composition due to high costs, confirmed that the observed proportions of gap sizes followed the pattern that we predicted using the virtual rock pile method. These results indicate that computer modeling can be used to predict the true proportions of gaps; however, the real proportions may be more variable than predicted. The computer modeling method also appears to have a tendency toward underestimating the number of small and optimal gaps, and overestimating the number of large gaps. The difference in error between the large composition and the mixed/round compositions indicates that predictions may be more accurate for smaller rocks than for larger rocks.
The use of rock models produces a more accurate model than would be obtained using simpler shapes such as cubes. However, this also makes the simulation computationally intensive as it means modeling a large number of polygons and their interactions. A limitation of Unity is that it cannot easily model concave rigidbodies. All the rock models were therefore convex, although real rocks often have some concavities. Despite this, our predictions approximately matched the gaps present in real CT-scanned rock piles (with rocks containing concavities). The accuracy of predictions for the entirely convex river rocks and the mixed compositions was fairly similar, suggesting that modeling the rocks as convex did not have a substantial influence on the accuracy of the prediction.
As it is impossible to delineate discrete gaps within the network of interstitial space, and as a contiguous space will vary in its dimensions depending on the point at which it is measured, we took an approach that estimated overall proportions of gap space of different sizes rather than attempting to count individual spaces. Our approach, however, does not directly provide information about the connectivity of gaps. Even a large gap could be completely cut off and inaccessible, for instance at the entrance to the pile. Rock-dwelling animals are often highly flexible and can bend around corners within piles easily; modeling this movement is more challenging than measuring gap sizes. There is therefore an underlying assumption to this model that the proportion of theoretically accessible space corresponds to the proportion of actually accessible space.
Using predictions from the model to design refuges Our method of protecting prey from predators relies on exploiting differences in their size. It is thus important to consider how different sizes of skinks will be affected, including both pregnant and large adults, as predation rates may be affected by size. Specifically, some evidence suggests that larger skinks are more vulnerable to rodent predation (Towns 1991, Newman 1994, Nelson et al. 2016, which has been attributed to the idea that smaller lizards are better able to avoid predation in refuges inaccessible to rodents (Towns 1991). If this is correct, then creating rock piles with optimized interstitial spaces may mean a greater number of spaces sized for larger adults to evade predators, reducing the impact on larger adults. It also means that this approach may be particularly useful for increasing skink recruitment, as juveniles are able to make use of small interstitial spaces to evade predators. Other evidence indicates high numbers of mice are correlated with suppressed skink recruitment, suggesting that when neonates/juveniles are abundant and vulnerable, development of a search image by mice is easier and more rewarding (Wedding 2007). The presence of interstitial refuges may make it harder for mice to prey on neonates and reduce their likelihood of forming a search image during periods of neonate abundance.
Pregnant skinks, being larger, likely need larger retreats than nonpregnant skinks and are more at risk of predation due to the physical burden of clutch weight (Shine 1980). However, viviparous species (as all but one New Zealand species are; Cree and Hare 2016) may have an advantage as their extra mass is more malleable than the rigid eggs of oviparous species (Schwarzkopf et al. 2010). Pregnant skinks were not tested in the access experiment (Appendix S1). If the optimized interstitial spaces are not large enough to protect pregnant skinks, this could make them more vulnerable to predation and reduce recruitment. Recruitment may also be affected by the size of individuals within the population; larger females produce larger clutches (Cree and Hare 2016). These factors need to be taken into consideration when designing rock piles, and for this reason, it may be prudent to ensure the presence of some larger gaps within the pile rather than seeking to include only optimally sized gaps.
To protect skinks from mice predation, we posit that maximizing the proportion of large spaces and minimizing the amount of small spaces (while also maximizing optimal spaces) are the best approach, as a rock pile that is largely inaccessible to skinks is much less useful than one that is largely accessible to skinks with some areas that are inaccessible to mice. Erring on the side of providing habitat for pregnant and large skinks is also advisable because survival of pregnant skinks is necessary for population survival, and an inadvertent increase in predation pressure on large individuals may have unintended evolutionary consequences. Some Oligosoma species are known to be aggressive and territorial (van Winkel et al. 2018), so a large number of optimal gaps are necessary to sustain a large population. Although large gaps may still form a sizable proportion of the rock pile, having more optimal gaps reduces the size of the large gaps compared with piles of larger rocks, so there will be fewer spaces large enough for mice to build nests (i.e., an 11-mm gap is preferable to a 60-mm gap).
Within the modeled rock piles, as optimal space increases, small space also increases at approximately the same rate. This makes it difficult to maximize optimal gaps while minimizing small gaps, as all the mixes are compromises with no clear best ratio. However, we hypothesize that a rock pile composed of entirely larger rocks (20-40 mm) will be best as it has the lowest proportion of small gaps and is therefore likely to be the most accessible to skinks, and still has a reasonable proportion of optimal spaces (41%; Fig. 3).
As our findings indicate a narrow range of gap sizes that benefit skinks, an alternative approach could be to create an artificial structure with precisely sized gaps. A benefit of the rock pile approach we recommend is that there are a diversity and complexity of gaps which can encourage diversity of invertebrates (a food source) by providing variable habitat (Tokeshi and Arakaki 2012). Additionally, artificial rocks which mimic natural rock attributes (cavity geometry, thermal properties) appear to be preferable retreats for reptiles and invertebrates than simple concrete pavers Shine 2000, Croak et al. 2010). However, the alternative approach allows for precise control of the size of gaps, which could be an improvement over juggling ratios in rock piles. These two approaches could be compared experimentally.

CONCLUSIONS
We promote the use of 20-40 mm rocks in the creation of rock piles for conservation of skinks of a similar size to or smaller than Oligosoma aeneum and Oligosoma polychroma. We acknowledge this may not be a sufficient recipe to protect larger species of skink in New Zealand. This is a hypothesized best design and should be tested in experiments with skinks and mice.
Our computer modeling technique allows for the characterization and measurement of interstitial spaces that cannot be measured by other means, and can be generalized to systems beyond terrestrial rock piles. In aquatic environments, habitat complexity (including factors such as size and makeup of particulate substrate and arrangements of fractal habitat elements) is an important influence on ecological communities (Tokeshi and Arakaki 2012). These factors influence interstitial space availability, which is currently measured using physical tools which can be limited in their scope and v www.esajournals.org highly labor-intensive to use (McCormick 1994, Tokeshi and Arakaki 2012, Rogers et al. 2017, or as mean volume which is less informative (e.g., Levine et al. 2017). These types of studies could benefit from the virtual technique we designed here, which would enable easier, more informative characterization of interstitial spaces, and contribute to understanding of habitat complexity.
The techniques outlined here could be used to inform management of many species relying on refuge spaces for some aspect of their life history. In the current study, the skinks used were relatively small, allowing the size difference between skinks and mice to be exploited. This approach can be generalized to systems where there is a difference between two size classes (e.g., predator and prey) that can be exploited. Examples of ecologically important differences in accessibility can be found in predator-prey systems (e.g., smaller prey able to escape larger predators ;Towns 1996); in cases of sexual dimorphism (e.g., smaller female fisher [Pekania pennanti] able to access natal dens larger males cannot; Green et al. 2019); or where there are differences between juveniles and adults (e.g., smaller/younger individuals have different burrow preferences than older/larger individuals; Milne and Bull 2000).