Conducting Source-Sink Analysis

Overview

Landscape source-sink properties emerge from HexSim simulations without the need for assumptions about the spatial distributions of vital rates or of movement frequencies. Performing HexSim source-sink analyses does not necessitate that modifications be made to the landscapes influencing a simulated population. Instead, one or more map tessellations are used to sample movement or vital rates, but this sampling can be conducted without impacting the population.

The key steps involved in HexSim source-sink analyses are listed below:

Construct the simulation model.
Develop one or more patch maps.
Add machinery to record locations.
Run the simulation.
Generate HexSim reports.
Quantify source-sink behavior.

Steps 2, 3, 5, and 6 are described in detail below. The generation of patch maps, and the analysis of the Productivity and Projection matrix reports, will be made easier through the use of the software utilities available at this web site.

References to software available in the Utilities section are highlighted.

Building Patch Maps

The term patch map is used here to represent any Hexmap that is comprised of hexagons having integer-valued scores. A patch, in this context, is simply the collection of all hexagons having the same score. Thus, a patch may consist of strictly adjacent hexagons, or it may be partly or fully disconnected.

Source-sink analysis can be performed with one or more patch maps. When multiple patch maps are used in the same simulation, and if they are all sampled at the same point in the simulation, then each map will provide a unique but equivalent assessment of the source-sink structure.

Any collection of one or more patch maps may be used to conduct a HexSim source-sink analysis. Each patch map must contain at least two patches. The source-sink analysis described here assumes that each patch map is static. That is, the patch maps used in these analyses must not change with time as a simulation runs.

Linking Patch Maps to Models

Source-sink dynamics can be captured if each individual records which patch it is in, and if individuals update this information each time they move. An example patch map is shown to the left. The labels were added outside of HexSim.

This image includes 9 habitat patches embedded within a non-habitat matrix (shown in black). The modeled population was excluded from the purple areas in the image. For this patch map to be useful, the simulated population will require an accumulated trait with 10 values. The black portions of the map were scored zero, and the colored patches were scored as labeled. Thus the 10 trait values shown in the figure below correspond to the matrix region plus the 9 habitat patches.

The Accumulated Trait data shown to the right illustrates a location trait actually used to capture the location information for individuals occupying this patch map. This simulation model had a single Movement event. An Accumulate event was placed immediately after this Movement event. The Accumulate event used an Individual Locations updater function (linked to the patch map and the location trait) to set the location trait value equal to the patch each individual occupied. Every individual had to be in either the matrix, or one of the 9 patches.

Users who want to conduct source-sink analysis using multiple sampling schemes simultaneously need only replicate the procedures described here for each separate patch map. And each patch map used must be paired with its own dedicated accumulated trait. A separate Individual Locations updater function must be added for each patch map. However, all of these updater functions may be housed within the same Accumulate event.

Example scenarios illustrating these steps can be found in the sample workspaces that are available from this website.

The Reports

HexSim Reports and Tallies are described in detail in the User's Guide. Here, we will discuss the Productivity and Projection Matrix reports because they are used in HexSim source-sink analyses. Refer to the User's Guide for a more thorough presentation of these tools.

The Productivity report can be used to build a trait-stratified summation of birth and death records. The report uses information from the log file to create a CSV file with four columns: trait index, births, deaths, and births - deaths. If the report has been stratified by location trait, then the data in the births - deaths column can be used as a measure of source sink quality. An example for the above patch map is shown to the left.

This data was gathered after the simulation had reached steady-state. When a system is at steady-state, it can probably be assumed that, for any given patch, births - deaths is approximately equal to emigration - immigration. The Productivity report provides this birth - death information, and HexSim's Projection Matrix report can be used to obtain values for emigration - immigration. Users decide which metric they will use to quantify landscape source-sink structure.

The Projection Matrix report shown to the right can be used to create a matrix capturing the likelihood that individuals, having one combination of trait values in time step N, will transition to another trait value combination in time step N+1. If a single trait is used to stratify this report, and if that trait contains stage class information, then the report will produce a Leslie matrix (sort of).

If a single location trait is used to stratify a Projection Matrix report, then the result will be a matrix showing the probabilities or frequencies of moving between any two landscape patches. For the analysis here, users must construct location-stratified Projection Matrix reports containing count data. This is done by making sure the "Output Raw Count Data" check-box has been selected when the report is produced. The above example Projection Matrix report was produced this way, so it corresponds to the above patch map. Both the rows and columns in this report correspond to the ten patches, using the order found in the corresponding location trait. The column indicates the starting location, and the row indicates the ending location. Thus, the value of 1037 in cell B-4 above indicates that 1037 individuals were observed moving from patch 1 to patch 3 (the first patch is labeled zero).

For location-stratified Projection Matrix reports, each column sum provides the total number of emigration events out of a patch, and each row sum provides the total number of immigration events into a patch. Thus, our second estimate of source-sink value can be obtained from the Projection Matrix report as a difference of row and column sums. The matrix shown above was also gathered after the simulation had reached steady-state.

The image to the left provides a comparison of these two measures of source-sink quality for this system. For each of the 10 patches, the birth - death (B - D) values were obtained directly from the Productivity report. The emigration - immigration (E - I) values were obtained from the Projection Matrix report. A plot of columns B vs. C produces a very straight line.

The final step in the process is typically to produce a Hexmap from the source-sink data obtained using the Productivity or Projection Matrix reports. Conceptually, this involves a simple mapping of the desired source-sink metric's values back onto each hexagon in the patch map. For example, patch 7 (see the Hexmap shown above) is composed of 59 individual hexagons that are all scored 7. If the source-sink value obtained from E-I was to be assigned to this patch, then all 59 of those hexagons would have to be changed to have a score of 3916, since that is the source-sink value for patch 7 seen in the right-most column of the above table. In practice, this can be done most easily using a small computer program that reads a Hexmap (exported to a CSV file) and either a Productivity or Projection Matrix report. Two such computer programs are available in the Utilities section of this web site.

An Example

The example presented below is an extension of the work described in Schumaker et. al. 2014. Please refer to that paper for details on the northern spotted owl simulation model referenced here. The Schumaker et. al. paper performed a source-sink analysis using two spatial stratification schemes, both of which were very course. This example illustrates how the methodology described above can be exploited to obtain a highly detailed map of source and sink quality distributed across the range of the northern spotted owl.

The concepts and analysis described in this example were developed by Nathan Schumaker. The first step in the process involved the generation of a suite of 25 different patch maps.

For this purpose, I developed a utility called build_hexmap_hexagons.c that constructed regular Hexmap tessellations using hexagonal-shaped patches. Each patch was comprised of multiple "atomic" hexagons. Patch area varied from 397 to 7651 hexagons, and the corresponding number of patches ranged from 803 to 58 (see the figure above-right).

Each of these 25 Hexmaps was originally produced as a CSV file. These maps were then clipped to the outline of northern spotted owl's range, using the utility called clip_hexmap.c, and then re-indexed using the renumber_patches.c utility. All hexagons outside the range of the owl were assigned a score of zero (shown in black). This clipping process caused some patches to have less area than others, since some atomic hexagons along the range edges were set to zero. A close-up image of the patch map structure is provided in the image to the left.

For every patch map, a corresponding accumulated trait was added to the spotted owl simulation model. And all of these traits were updated at the same spot in the event sequence to ensure that each map produced an equivalent quantification of the overall source-sink dynamics generated by the simulation. Once these steps were completed, the intent was to run 100 replicate simulations, and then to generate a single Projection Matrix report for each of the 25 patch maps. The reports were to be constructed from a combined log file that contained the results from all 100 replicates.

A big problem arose at this point in the process. The northern spotted owl simulation model depended on a suite of life history traits that together produced a total of 864 trait combinations. This work was being done on a 64-bit computer, so the maximum allowed number of trait combinations was 2^64, which roughly equals 1.85 x 10^19. This is a large number, and the initial 864 trait combinations were nowhere near the limit. However, when the accumulated traits corresponding to all 25 patch maps were added to the simulation, the result was a model with 1.00 x 10^67 trait value combinations. This massive number was 45 orders of magnitude too large for the 64-bit computer to handle.

The solution was to break the simulation up into four separate models that each used a subset of the 25 patch maps. Doing so, it was possible to keep each model's total number of trait combinations below 2^64, while collectively making use of all 25 patch maps. It was then necessary to ensure that each of the four models produced 100 identical replicate simulations, with the only differences between them being the location information used to generate the Projection Matrix reports. Thus a collection of 100 random number seeds was assembled to ensure that each of the four groups of 100 replicate simulations were identical. Random number seeds can be supplied to HexSim when the model engine is called from the Windows Command Line tool. The assignment of individuals to patches does not result in any calls being made to the random number generator, so truly identical replicates could be obtained this way using different patch maps (this assertion was thoroughly tested).

Each of the 4 models was used to run 100 replicate simulations.

Each simulation built a log file, and each collection of 100 log files was combined into a single huge file.

The 4 combined log files were each 0.25 Terra-bytes in size. These 4 files were used to generate 25 Projection Matrix reports, representing the 25 patch maps.

Next, the 25 Projection Matrix reports were used to create 25 different source-sink maps. This was accomplished using map_projection_counts.c.

Each source-sink map was then re-scaled so that the sink hexagon scores ranged from -100 to 0, and the source hexagon scores ranged from 0 to +100. 24 of the 25 source-sink maps created this way are shown in the figure to the right.

The color scale used to display the images was altered to create images with higher contrast. This color scale is available from the Accessories page of this web site (see Red-Yellow.pal) The image built up from the largest patches was left out to keep the above image compact. Neutral areas (those that are neither sources nor sinks) are displayed in black in the above figure. Much of the Pacific Northwest is black in the image because, at steady-state, the simulated owl populations had dropped out of this portion of the landscape.

Finally, the 25 uniformly re-scaled source-sink maps were summed to produce a single multi-scale model of northern spotted owl source-sink dynamics across the species' range (see image). A comparison of the combined multi-scale image (above-right) with the individual single-scale source-sink maps (above-left) suggests that no single scale of analysis may be sufficient to produce a defensible source-sink assessment.

Below is a summary of the steps used to obtain the multi-scale composite source-sink map depicted above:

Develop a suite of patch maps representing a spectrum of different spatial scales. Ideally, trim and renumber these maps so that no patches fall completely outside the study area (e.g. species range).
Add an accumulated trait to the simulation for each patch map. The number of trait values will be equal to the number of patches, or number of patches + 1 if the matrix outside the patches is to be considered. Be sure the total number of trait values in your scenario does not exceed 264 (on a 64-bit computer) or 232 (on a 32-bit computer).
Use Individual Locations updater functions (within Accumulate events) to set individual's trait values based on the patch maps. This will allow you to stratify Productivity or Projection Matrix reports based on the patch map-specific traits.
Run the simulation, making sure to obtain a log file.
Construct a separate Productivity or Projection Matrix report for each patch map used.
Use the map_productivity_report or map_projection_counts utilities to convert the reports into Hexmaps. These programs are available from the HexSim website under Resources / Utilities. The utility will need to be run once for each patch map, using a copy of that patch map that's been exported as a CSV file. These utilities will generate a new CSV file that represents a Hexmap.
Re-scale each individual source-sink Hexmap produced in the previous step, so that its values fall between [-100 , 100]. The actual scaling values (e.g. +/- 100) aren't important, as long as every Hexmap produced in the previous step is mapped onto the same fixed range. Recall that the individual Hexmaps are assembled using the patch maps, with every hexagon comprising a patch being assigned the same value. The upper-limit on the value assigned to the patches increases in magnitude as patch size increases. In this step, we are ensuring that each single-scale analysis contributes equally to the final multi-scale map, regardless of patch size.
Add all of the re-scaled Hexmaps together. When doing this, be sure that each Hexmap has exactly the same number of rows, and that they are sorted the same way (i.e. from lowest to highest hexagon ID). Only add the score columns together, so that the Hexagon IDs are unchanged. The resulting CSV file should have a 1-line header, and two data columns (Hexagon ID and Hexagon Score).
Use HexSim's Display Spatial Data tool (found in the HexSim menu) to read the resulting multi-scale composite source-sink Hexmap. You may first want to edit the CSV file and clip the minimum and maximum values to enhance the image contrast. For example, if most source hexagons had a value falling between 1 and 100, but a few hexagons were scored 1000, then most source hexagons will be assigned similar colors. If the few hexagons scored greater than 100 were manually assigned scores of 100, then the Display Spatial Data tool will produce a more useful image. The same is true for sink hexagons, which have scores less than zero.
Optionally, export the multi-scale composite source-sink Hexmap to a bitmap. Then, consider swapping the color map with the one used to generate the image shown above. That color map is available from the HexSim website under Resources / Accessories. The file in question, named Red-Yellow.pal, can be used with XnView to modify the bitmap's color scheme.

This 10-step process was applied to a simulation model of Brazilian Maned Wolves to produce the source-sink image shown above. In this case, four patch maps were used to capture source-sink information at four different scales. This source-sink information was then quantified using four separate Productivity reports. Extreme source and sink values were clipped to enhance image contrast. Also, source values below 10 were set equal to 10, so that the source-sink CSV file could be concatenated with a mask file having values of 1 that provided the geographic context visible in the image (the gray study area). This allowed all neutral study area pixels (neither sources or sinks) to be assigned a single color.