Supplementary Material
Appendix I. Details of the Cellular Automata Model.
Although prior studies provide clear evidence for gross increases or decreases in the vital rates over spatial gradients, the precise functional forms of most parameters are not known. Therefore, we assumed the simplest functional relationships supported by available data, and then use the model to generate qualitative trends over an idealized two-dimensional landscape, rather than make specific numerical predictions. In the Discussion we consider the effects of choosing different functional forms for key parameters.
In the cellular automata (CA) model, a mussel bed is represented as a rectangular lattice of cells, where each cell is a potential site for the location of a mussel. Each cell is identified by its coordinate in a lattice, where the horizontal dimension represents an alongshore gradient in wave energy and the vertical dimension represents a gradient of tidal height. To facilitate comparisons with lattices of different sizes, the and coordinates are scaled to vary from zero to one. Each cell can take on integer state values representing mussel size. The state of cell in position at time is , where 0 is an empty cell, 1 is a cell with a newly recruited mussel, and is the asymptotic maximum size of a mussel at a particular coordinate, which varies with shore level and wave energy. Model dynamics are described as stochastic transitions among states, where the transition probabilities for any given cell are determined by the current state of the cell, the cell’s position in the lattice (gradient effects), and the states of the surrounding cells (“neighborhood effects”). The transition probabilities represent position-specific mussel recruitment, growth to larger sizes, size dependent predation, and background mortality, which includes mortality for physical stress.
Spatial patterns of mussel recruitment are perhaps the least understood aspect of mussel ecology. It is known, however, that M. californianus recruitment probabilities often increase with wave energy along a given shore level, reaching a maximum on wave beaten shores (Menge 1992; Robles 1997; Moya 2005). To estimate a 2-D surface of mussel recruitment, vertical transects sampling naturally occurring M. californianus recruits (1-5 mm long) were placed at regular alongshore intervals on a shore with a horizontal gradient of wave energy. Samples were made late winter of three successive years, the season with the highest recruitment of M. californianus on the Catalina site (Robles 1997). The recruitment surface showed a peak occurring in the mid-intertidal zone at the wave-exposed extreme and decreasing towards lower wave energies and more extreme shore levels (Figure S1). The CA model assumes that recruitment rates to empty cells peak at mid-shore levels and high wave energies. Along the vertical axis, recruitment rate declines from the peak mid-shore value to higher and lower shore levels according to a Gaussian. Along the horizontal axis, recruitment rate declines exponentially with decreasing wave energy. The gradient effects in recruitment were summarized by the following function:
, (A1)
where , and are minimum and maximum recruitment rates along the wave energy gradient, is the shore level where recruitment rates reach a peak, and is a parameter which specifies the rate of recruitment decline with shore level as one moves away from the mid-tidal peak.
Mussels have indeterminate, environmentally influenced growth (Sebens 1987). Growth rates and maximum (“terminal”) sizes increase with longer immersion times (i.e. lower shores levels; Dehnel 1956; Garza 2005), and higher nutrient flux (i.e. higher wave energies; Leigh et al. 1987; Dahlhoff and Menge 1996). In the model, mussel growth rates were expressed with the von Bertalanffy function:
. (A2)
The expected growth increment of a mussel, , is proportional to the difference between the size of a mussel, , and its maximum asymptotic size, , where is the growth rate. Maximum asymptotic size depends on flow rate and immersion time and, therefore, increases with increasing wave energy and decreasing shore level. The CA model uses a simple product of linear trends along the horizontal and vertical dimensions:
, (A3)
where and are the minimum and maximum asymptotic sizes on the lattice.
The CA model includes two possible sources of mussel mortality: background mortality and predation. For simplicity, the same low rate of background mortality is assumed for all mussels independent of their size and location in the lattice. Mortality rates due to predation are the product of the per capita attack rate of predators, , and the predator density, . Predation rates decrease with lengthening periods of tidal emersion (i.e. vertically towards higher shore levels; Paine 1974; Menge 1992; Robles et al. 1995, Garza 2005) and greater hydrodynamic stress (i.e. horizontally from sheltered to wave-exposed areas; Menge 1983). The CA model uses a simple linear trend from a maximum predation rate at the lowest wave energy and lowest shore level to a minimum rate at the highest wave energy and highest shore level:
. (A4)
All other rates are assumed to be constant within the lattice.
Rates of recruitment and predatory mortality are modified by neighborhood effects. Local densities of mean mussel biomass were computed for neighborhoods of size as a simple average of the size of the mussels in the grid of cells centered at :
. (A5)
At the sides and corners of the lattice, the local densities were calculated as the mean of the subset of cells within the neighborhood.
Our observations and field experiments (Moya 2005) indicate that while M. californianus can settle out on any rough surface, and recruitment rates are higher under adults than algae or bare rock (Robles unpublished data). The CA model assumes that rates of recruitment to empty cells increase exponentially with the density of surrounding mussel biomass:
, (A6)
where is a parameter representing the strength of the recruitment neighborhood effect, is the size of the recruitment neighborhood, and is given by equation (A1).
Experimental studies (Bertness and Grosholtz 1985; Robles et al. 2009; Fong and Robles, unpublished data) indicate that younger, smaller mussels are shielded from predators by the larger less vulnerable members of an aggregation. Thus, the likelihood of predation decreases as a mussel becomes surrounded by larger, less vulnerable individuals. The CA model assumes that large mussels are more resistant to predation and that mussels surrounded by large conspecifics are less susceptible to attack. This neighborhood effect is modeled as a negative exponential decrease in predator attack rate with an increase in the density of neighborhood mussel biomass. The total mortality rate of a mussel in cell is given by
, (A7)
where is the strength of the predation neighborhood effect, is the size of the predation neighborhood, is the background mortality rate, and is given by equation (A4).
Experimental evidence shows that sea stars congregate in episodes massive recruitment of small mussels and disperse as these preferred prey are reduced in abundance relative to the larger mussels (Robles et al. 1995). This aspect of predation constitutes a numerical response. In the CA model, predator density is an additional global variable. The density of predators is modeled as an immigration-emigration process. It is assumed that predators enter the system at a constant rate and exit the system at a per capita rate that is inversely proportional to the overall per capita consumption rate of prey:
, (A8)
where the sum in the denominator is over all the cells in the lattice. The parameter is the constant of inverse proportionality between mussel biomass consumed and rate of emigration.
In the CA model, mussel recruitment, growth, and mortality are treated as stochastic events. For mussel recruitment and mortality, this was done by choosing an iteration time step and assuming rates are constant within that time step. This assumption yields an exponential function for the probabilities. The probability of a recruitment event into an empty cell is given by
. (A9)
Similarly, the probability that a mussel in cell dies is given by
. (A10)
Each mussel that dies is assigned a cause of death: background mortality with probability or predatory mortality with probability . The total biomass of mussels consumed by predators in each time step, , is used in the computation of predator emigration rates (see below). If a mussel escapes a random mortality event, then it may experience a growth event in that time step. While this assumed ordering of survival followed growth seems arbitrary, from a practical consideration, growth occurs slowly on the time scale of our model iterations, so that this assumption is of little importance. A truncated Poisson distribution is used to grow mussels:
(A11)
where is the expected growth increment over the time interval .
Predator immigration and emigration are handled as a stochastic birth-death process. However, predator movements occur on a much faster time scale than changes in mussel biomass. We, therefore, found it necessary to implement a second, shorter time step of for predator immigration and emigration events. During the time interval , predators randomly enter and leave the system every based on mussel biomass densities from the previous time step; i.e. mussel densities are not updated during the time steps of size . Without this assumption, the expected number of predator immigrants is , which may, in fact, exceed the theoretical equilibrium number of predators unless is very small. On the fast time scale, the predator immigration probabilities are given by a Poisson distribution:
, (A12)
where is the number of predators entering the system. For emigration, the probability of an individual predator leaving the system was calculated using equation (A8) where the summation in the denominator was replaced with , the actual biomass of mussels eaten during the previous time interval . The probability that a predator emigrates during the time interval equals
, (A13)
and the probabilities for the number of emigrating predators is given by the binomial distribution:
, (A14)
where is the number of predators leaving the system and is the number of predators at the beginning of the time interval (product of the lattice size and the predator density). At the end of each time interval , the new predator density equals
. (A15)
Equation (A15) is iterated times for every time step.
Parameter values used in the simulations appear in Table S1.
Table S1. Variables and Parameters of the Cellular Automata Model.
Time
/ days (d) / Time
/ 1 – 8 d / mussel bed time step; varies†
/ d / predator time step; varies†
/ (unitless) / number of predator time steps per mussel bed time step; varies†
Space
area / 1 cell / each cell holds a single mussel
/ (unitless) / horizontal location of a mussel in the lattice, scaled 0 to 1
/ (unitless) / vertical location of a mussel in the lattice, scaled 0 to 1
/ 10000 or 16000 / width (number of cell widths) of the mussel bed lattice; 16000 cells on single within-site lattice; 10000 cells on each between-site lattice
/ 250 / height (number of cell heights) of the mussel bed lattice
Mussel Size
/ mm / size of a mussel in position
/ 1 (unitless) / neighborhood radius; number of cells in each direction
/ mm / mean size of mussels within radius r of position (Eq A5)
Mussel Recruitment
/ d-1 / mussel recruitment rate into empty cell at position without neighborhood effects (Eq A1)
/ 0.0 d-1 / minimum recruitment rate along the wave energy gradient
/ 1.0 d-1 / maximum recruitment rate along the wave energy gradient
/ 0.5 (unitless) / tidal height location of mussel recruitment peak
/ 0.01 (unitless) / shore level width of mussel recruitment curve around peak
a / 0.05 mm-1 / enhancement of recruitment with neighboring mussels
b / 1 (unitless) / recruitment neighborhood radius
/ d-1 / mussel recruitment rate into empty cell at position with recruitment neighborhood effects (Eq A6)
Mussel Growth
/ mm / asymptotic size of a mussel in position (Eq A3)
/ 30 mm / minimum asymptotic size of a mussel within lattice gradients
/ 200 mm / maximum asymptotic size of a mussel within lattice gradients
/ 0.0007 d-1 / decrease in mussel growth rate with size
/ mm d-1 / mean size change per unit time of mussel in position (Eq A2)
/ mm / expected mussel growth increment
Mussel Mortality
/ 0.0001 d-1 / background per capita mussel mortality rate
/ d-1 / total per capita mortality rate of a mussel in position (Eq A7)
/ cell pred-1 d-1 / predator attack rate at position (Eq A4)
/ 0.0 cell pred-1 d-1 / minimum predator attack rate within lattice gradients
/ 1.0 cell pred-1 d-1 / maximum predator attack rate within lattice gradients
c / 0.04 mm-1 / decrease in predator attack rate with mean mussel size
d / 1 (unitless) / predation neighborhood radius
Predator Movement
/ pred cell-1 / predator density
/ 0.01 pred cell-1 d-1 / predator immigration rate
/ pred / number of predators immigrating into the lattice (Eq A12)
/ d-1 / per capita predator emigration rate (Eq A8)
/ 5.0 mm pred-1 d-2 / predator emigration constant per unit prey consumed
/ mm / total mussel biomass consumed during previous time interval
/ (unitless) / probability that predator emigrates during interval (Eq A13)
/ pred / number of predators emigrating from the lattice (Eq A14)
† Time steps were varied using an optimization procedure.