J. theor. Biol. (1992) 159, 111-133
Is the Crown-of-thorns Starfish Degrading the Great Barrier Reef? R. M. SEYMOUR~" AND R. H. BRADBURY~
t Department of Mathematics, University College, Gower Street, London WC1E 6BT, U.K., and ~ National Resource Information Centre, GPO Box 858, Canberra A C T 2601, Australia (Received on 7 November 1992, Accepted in revised form on 9 May 1992) The phenomenon of crown-of-thorns starfish outbreaks on the Great Barrier Reef (GBR) is treated on a large scale as analogous to a disease (starfish outbreaks) spreading through a community of susceptibles (individual reefs). A simple (epidemiological) model is found which well represents the extant data on starfish abundance for the central sector of the GBR. The form of this model suggests that some longterm degradation of reef community structure is taking place. The model is also used to predict the gross outbreak behaviour for the central sector of the GBR up to the year 2000. A comparison with the (sparse) data for the northern and southern sectors of the GBR is also made which casts some doubt on a hypothesis of Kenchington (1977) concerning the primary locus of outbreaks.
1.
Introduction
Infestations of coral reefs in the Indo-Pacific region by very large numbers of crownof-thorns starfish (Acanthaster planci) have been observed since the early 1960s, and have generated much controversy concerning their possible causes (Moran, 1986). In spite of a huge research effort spanning more than 20 years, the spatial and temporal extent of these outbreaks is poorly documented in general. Many of the reasons for this paucity of good data are discussed in Moran (1986). Documentation efforts have been pursued most consistently in the central sector of the Australian Great Barrier Reef [GBR, (Lats. 16°-20°S)], and reasonably good data do exist for this region (i.e. a reasonably extensive sequence of sample points of reef states, based on reasonably large sample sizes, exists--see section 2 below). A lesser amount of data also exist for the northern and southern sectors (Lats. 10°-15°S and 21°-24°S). These data apparently indicate southward moving waves of infestation, presumably driven by the East Australian current, via transportation of starfish larvae (Reichelt et al., 1990). In this paper this Acanthaster phenomenon on the GBR is considered as analogous to a disease spreading through a community of susceptibles of fixed size, in which the "susceptibles" are individual reefs, and the "disease" is starfish infestation (Frauenthal, 1980; Beeker, 1989). The data are found to be well represented by a simple (epidemiological) model with a constant recovery rate for reefs, but a timedependent infection rate. The latter incorporates a short-term periodicity (approximately 12 years), but perhaps more interestingly, a much longer term exponential III
0022-5193/92/2101 ! I +23 $08.00/0
© 1992 AcademicPress Limited
112
R. M.
SEYMOUR
AND
R.
H. BRADBURY
decay. Both these features are interpreted as representing coral processes, with the short-term periodicity providing clear evidence of a predator/prey type cycle, and the decay term indicating long-term degradation of coral community structure. The findings for the central sector of the GBR are also related to the (lesser amounts of) data which exist for the northern and southern sectors. This comparison indicates that a currently popular model of a wave of infestation generated in the northern sector and travelling southward (Reichelt et al., 1990), does not seem to be strongly supported by the extant data. Finally, the main model obtained from the central sector data is used to predict the large scale outbreak behaviour over this sector to the year 2000. The paper is arranged as follows. In section 2 the data set used is discussed, and confidence weightings based on Bayesian statistics are obtained. The theory used is given in the Appendix. In section 3 the main model is constructed and fitted to the central sector data (using two different weightings). These results are compared with those obtained from some possible alternative models in section 4. The comparison with the northern and southern sector data is made in section 5. Finally, discussion, interpretation and prediction are given in section 6.
2. The Data Set
The data set used was extracted from a database of historical records of Acanthaster abundance obtained from surveys of samples of reefs over the 24-year period from 1966 to 1989 (Bradbury & Mundy, 1989). Two state categories were used, 0 indicating low starfish abundance on a reef (uninfected), and 1 indicating starfish abundance of outbreak proportions (infected). Within each yearly sample of reefs, the number of reefs surveyed in each state were counted. Data of the above type exist, from some region of the GBR, for each of the 24 years. However, largely for logistical convenience, the GBR has been divided into three distinct zones; the northern sector, the central sector and the southern sector, and a breakdown of the data by sector shows heavy sampling bias towards the central sector (the full data sets for each sector are given in Tables l, 2 and 3). Thus, there are 22 out of a possible 24 sample points for the central sector, with 16 for the northern and 17 for the southern sectors. Further, more than 40% of the total number of reefs sampled in any 1 year come from the central sector in 20 out of the 24 years, while more than 60% come from this sector in 12 out of the 24 years. The full sampling distribution for the proportion of reefs coming from the central sector is given in Fig. 1. It is mainly due to this bias that the modelling effort is concentrated on the central sector. While it would no doubt be possible to eliminate additional sampling bias by considering narrower latitude bands, this process would fragment the data still further. Concentration on the central sector eliminates what is almost certainly the most significant bias while still retaining a useful amount of data. A further reason for dividing the data into sectors in this way will be considered in connection with the "Kenchington wave" hypothesis in section 5.
DEGRADATION
BY T H E C R O W N - O F - T H O R N S
113
STARFISH
TABLE 1
Data for the central sector of the GBR n E s CI C2
CL n E s CI C2
CL
0 45 30 0.86 0-97 0.79 0.53
1 I0 8 0-6 0-78 0.96 0.52
2 4 2 0.37 0.53 0'85 0.15
3 35 14 0.79 0-94 0-56 0.25
4 69 18 0.94 0.99 0-37 0-17
5 43 23 0-82 0.96 0.68 0.39
6 26 19 0-77 0.92 0.87 0.55
7 25 13 0-7 0-89 0.7 0.34
8 25 12 0.7 0'89 0.67 0-3
9 --0 0 ---
10 15 2 0.74 0.87 0.35 0'03
11 --0 0 ---
12 2 I 0.3 0.44 0.91 0-09
13 5 4 0-47 0-65 0.98 0.41
14 35 7 0.86 0.97 0.35 0.1
15 29 13 0-74 0.91 0.62 0.28
16 28 13 0.73 0.91 0-64 0-29
17 74 28 0-93 0.99 0.49 0.27
18 87 29 0-96 !-0 0.44 0.24
19 86 14 0.98 1-0 0-25 0.1
20 61 18 0.92 0-99 0.42 0-19
21 62 9 0.96 I-0 0-25 0.07
22 71 15 0-96 !.0 0.32 0.13
23 70 9 0.98 1.0 0.22 0'06
n = y e a r ( 0 = 1966); ~Z=number of reefs sampled; s = n u m b e r of reefs in sample which are in state 1 ; C I = confidence estimate with ~ = 0-1 ; C2 = confidence estimate with ~ = 0-15; C L = 95% confidence limits for s / Z (see text for explanations); - - = no data.
TABLE 2
Data for the northern sector of the GBR n E s C2
CL n g s C2
CL
0 21 6 0.89 0.49 0-13
1 12 5 0.75 0.68 0.19
2 --0 ---
3 --0 ---
4 16 4 0.85 0.48 0"09
12 --0 ---
13 5 4 0-65 0-98 0.41
14 50 13 0-98 0.39 0.15
15 9 8 0-78 0-99 0'6
16 22 7 0-89 0-52 0-16
5 . .
6 . .
0 . . 17 34 I1 0.95 0.49 0.19
. .
. .
. .
. .
0 . .
18 37 6 0-98 0.3 0.07
7
8
0
0
19 66 1 1.0 0.07 0
20 20 0 0-97 0.13 0
9 67 0 1.0 0.04 0
10 --0 ---
11 3 3 0.48 1.0 0-47
21 32 0 1-0 0.09 0
22 26 1 0-97 0.17 0
23 5 0 0-62 0.39 0
Designations as in Table 1.
It will be noticed that the sample sizes in Table 1 vary considerably from year to year (with no reefs sampled in years 9 and 11, up to 87 reefs sampled in year 18). Ifxn = Sn/Zn is the proportion of reefs sampled in year n which are in state 1 (infected), we would like to obtain some estimate of how accurately this reflects the actual proportion of central sector reefs which are in state 1 in year n. We believe that the most sensible way to approach this question is from a Bayesian point of view (Howson & Urbach, 1991). Thus, if an is this actual proportion, and we assume random sampling in the central sector together with uniform prior distributions for both an and sn, then we can ask: What is the probability that la,-xnl < 4, for some given error level 4? The details of the application of Bayes's theorem to obtain the relevant distribution are given in the Appendix, and the results for ~ = 0" 1 and 0.15 are given
114
R.
M. S E Y M O U R
AND
R. H . B R A D B U R Y
TABLE 3
Data for the southern sector of the GBR n £ s C2
CL
0 --0
I 3 0 0-48
2 --0
--
0-53
--
0
--
--
n E s C2
CL
12 14 0 0-91 0.18 0
13 9 3 0-71 0-63 0-11
14 1 0 0-28 0-78 0
3 14 4 0.81 0.53 0.11
4 3 1 0-51 0-77 0.04
15 3 1 0.51 0-77 0"04
16 9 4 0.68 0.73 0.18
5
6
. .
. . 0
. .
. .
. .
. .
0
. .
. .
17 31 7 0.95 0.39 0.11
18 34 5 0.97 0.29 0-06
7
8
0
0
19 40 6 0-98 0.28 0"06
20 28 3 0.96 0.25 0"03
9 91 17 1.0 0.27 0.12
10 2 0 0.4 0-63 0
11 I 0 0.28 0-78 0
21 33 4 0.97 0.26 0-04
22 --0 ---
23 3 0 0.48 0.53 0
Designations as in Table 1.
6 4
0
X 0"68
Iol r
0.13
I0
0 4
I
0"26
0'41
I
l
20
30
i 40
50
60
70
80
90
I00
] I loO l l 0'71
0.47
0"53
''
0.55
0"56
0.49
"1
0.9
FIG. l. Sampling distribution for the proportion (x) of reefs from the central sector (table); n = y e a r (0 = 1966); median = 0.62. The chart gives frequencies in each 10% range.
in Table 1. Notice that, with ~ =0.15, only years 2, 9, 11, 12, 13 fall below the 75% confidence level, so that 19 good data points are left if rejection is below this level. However, if this level of confidence with ~ = 0- l is wanted, nine data points would have to be rejected, while the same 19 data points as above can be retained if rejection is below the 55% level. So as not to be forced to reject too much data, or to have too low a level of confidence in much of what is retained, ~ = 0.15 was chosen in the sequel. Similar estimates for the northern and southern sector data are given in Tables 2 and 3. The procedure described above allows the introduction of various weightings for the data points, and in subsequent sections two such weightings will be used; a coarse weighting and a fine weighting. The coarse weighting is that described above: with =0-15, weight 0 is assigned to those points below the 75% confidence level, and
DEGRADATION
BY T H E C R O W N - O F - T H O R N S
STARFISH
115
weight 1 to those above it. The fine weighting is more subtle, and does not reject or accept data points so brutally. A one-sided normal curve is considered (from - c o to 0; total area= 1) with 5% tail at -0.25 (cr=0-25/1.96). With ~=0"15, any data point with confidence level <50% is assigned weight w0= tail area from - 0 . 5 to - c o . The remaining area under the normal curve is split into five parts, giving weights wk=area from [0-1 x ( k - 1 ) - 0 . 5 ] to [0-1 x k - 0 . 5 ] , 1 < k < 5 . A data point is then assigned weight Wk if its confidence level lies between [ 5 0 + 1 0 x ( k - 1 ) ] % and [50+ I0 x k]%. The resulting weighting is shown in Table 4. Clearly, many other weighting schemes are possible, and the effect of different choices will be discussed in section 6. TABLE 4
Fine weights C2
_<50%
>50% _<60%
>60% _<70%
>70% _<80%
>80% _<90%
>90% _<100%
w
0.0001
0"0016
0"0170
0.0982
0"3162
0"5669
C2=confidence level (%) with ~=0.15 (see text for explanation).
Finally, the sampling distribution described above may be used to obtain 95% confidence limits for the data points (see Appendix). These limits are given in Tables 1, 2 and 3. 3. The Main Model
The simple epidemiological model we consider is of the form
x,+t = x,-pn(1 ---,0). x. +pn(0 ---, 1). (1 - x , )
(1)
where x, is the fraction of reefs sampled which were in state 1 in year n, I is a latency period representing the natural time scale over which outbreaks develop, p,(1 --, 0) is the fraction of reefs in state 1 which revert to state 0 (recover) between years n and n+l, and p,(O~ 1) is the fraction of reefs in state 0 which become infested during this period. As it stands, (1) is merely tautological, and only acquires content when substantive assumptions are made concerning l,p,(1 ---,0) andp,(0 --* 1). In this and the following section several such assumptions are considered. Firstly, according to growth studies, A. planci reaches average adult size (~300 mm diameter) towards the end of its second year of life (Moran, 1986; Birkeland, 1989). In view of this, it is reasonable to take the latency period 1= 3 years, this being the time span (to the nearest whole year) from the settlement of larvae to the achievement of average adult size, which is the first stage at which outbreaks have usually been observed (Moran, 1986). Our second key assumption is that p,(l --, 0) is governed by features of Acanthaster life history; for example the ageing and death of an infesting starfish cohort. Field observations indicate that outbreaks form and disperse relatively rapidly (lasting a
116
R.M.
S E Y M O U R A N D R. H, B R A D B U R Y
few y e a r s at m o s t ) , a n d are usually f o u n d to be u n i m o d a l with respect to size class
(and therefore, presumably, age class--Moran, 1986). For this reason, we make the simplest possible hypothesis and take p,(1 ~ 0 ) = q, a constant. However, in spite of the success achieved with this assumption, it must be treated with some caution. Writing c = 1 - q , (1) may be transformed to
z~(c) = (x.+t-cxn)/(1 -x~) =p~(0 ~ I)
(2)
where both c and p.(O~ 1) are required to be in the range 0-1 [by definition--see (1)]. With 1= 3, and using the coarse weighting described in section 2, Table 1 gives 13 pairs (x., x.+/) for which good data exist. Evaluating zn(c) using this data shows that there is no allowable value of c for which all these values lie between 0 and 1 ; for small c and low values of n, z.(c) is >I, while for larger values of c, z.(c) is always <1, but some values then become negative. Thus, there may be significant variation in the recovery rate p.(1 ~ 0) over time. This point will be taken up again in sections 4 and 6. In contrast to the recovery rate, we assumed that the infestation rate p,,(O--, 1) is largely influenced by underlying coral processes occurring in response to starfish outbreaks, and which determine how susceptible a particular reef is to host an outbreak. This effect could not be modelled directly due to lack of suitable coral data, and thereforep. (0 --* 1) was allowed to be time dependent. To obtain an impression of the gross behaviour of this function, a regression line is fitted to the points z.(c) in (2)--evaluated from the data of Table l--using least squares estimation; i.e. by mlmmlzmg $2 = E w . [ z . ( c ) - a n - b] ~ n,
where w~ is the weight for year n. Thus, the following is obtained
a=Ac+B
(3)
where
A = {(E w.)(E w . n a . ) - ( E
w.n)(E woa.)}/A
B = {(~ w.n)(~, w.b.)- (~. w.)(~, w.nb.)}/A A = {(E w~)fE wnn2) -- (E w~'0 2}
a.=x.+t/(1 - x . ) ;
b.=x./(l - x . ) .
Using the coarse weighting (secton 2), this gives A =0.0849 and B=-0.0485. Thus, from (3), the regression has negative slope a for c<0'5713, and positive slope for c> 0.5713. However, if c is restricted to those values for which z.(c) is always nonnegative, then it must be that c~<0-26, and in this range, a in (3) is always negative. This shows that the overall trend in p.(0 ~ 1) should be to decline as n increases. In view of this, the preliminary guess for the form of this function is pn(0 ~ 1) = E . exp (-bn)
(4)
DEGRADATION
BY T H E
CROWN-OF-THORNS
STARFISH
117
To see how well this model captures the overall trend in the data, we proceed as follows. Using the coarse weighting, a system of running means from the data is formed
mk=(i~|O Wk+nXk+,)/(i~=|oWk+n)
(5)
(k >_0). With q= 9 this gives 16 points. For every choice of parameters c, E, b (with l= 3), and initial conditions 20, 21,22, the model (1) generates a complete sequence of points {2,10_
rhk(C,E, b, 20, 21,22) =
(Yo)/ 2k+~
q
(6)
(all weights= 1 in this case). The model means (6) may then be fitted to the data means (5) using least squares estimation for the six parameters c, E, b, 20,21,22. The value q = 9 was chosen as being large enough to smooth out short-term fluctuations in the data, but small enough to leave sufficiently many data points (residual variance estimated on 1 6 - 6 - 1 = 9 degrees of freedom). The result is shown in Fig. 2. The values c = 0- 1121 and b = 0.0774 confirm the long-term decline postulated in (4), and also that c is in the range for which z~(c) is always non-negative. However, this is not the whole story. With the values of c, E, b obtained above, the points E -~ . exp (bn). z,(c) may be plotted, and doing so indicates fairly clearly that there is also a periodic movement in the data. In view of this, a more refined guess for the form o f p , ( 0 ~ 1) is p~(0 ~ 1 ) = ( 1 / 2 ) . E . exp
(-bn).
{(1 + a ) + (1 - a ) . cos
[21r(n+k)/T]}.
(7)
The periodic term, ( 1 / 2 ) { - } , has period T, year-lag k, and varies between a and 1 (0 < a < 1) ; E is an overall amplitude scaling factor. Note that this form is as simple as we could expect for such a periodic term. With the extended form (7), every choice of parameters c, E, b, a, k, T (with l = 3), and initial conditions Xo, x~, x2, generates [via (1)] a time series {2,10_
S2= y w.(x.- 2.) 2. tl
(With the coarse weighting, the residual variance is estimated on 19 - 1 - 9 = 9 degrees of freedom.) The result for the coarse weighting is: c = 0.129, E = 1.2308, b = 0.0608, a=0.197, k = - 2 . 4 6 3 2 , T=I1.9938, 20=0.6771, 21=0"7908, 22=0"6255 (R z= 0- 9787; R 2 = 0- 8249, when corrected for degrees of freedom). This fitting procedure may also be performed using the fine weighting (section 2), and the result is shown in Fig. 3. Notice that all parameters have values commensurate with those obtained from the coarse weighting, but that the fit, as measured by R 2 or/~2, is not quite as good. This latter outcome is not surprising in view of
I I
I 2
I 5
I 4
I 5
I 6
I 7
I 8
I 9
1,, I0
I II
I 12
I 13
I 14
I 15
I 16
I 17
I 18
I 19
I 20
I 21
I 22
2:.'5
FIG. 2. Fit of running means for model to those of data with q = 9 (see text for explanation), c = 0. I!21, E=0.5698, b = 0-0774, ~?0=0.7158, :¢~ =0-9658, ~72=0.3705 [R2=0.9987; ~2=0.9897, when corrected for degrees of freedom]. • = d a t a point; 0 - - - - - 0 = model plot.
0'5
1.0
C
m
Z
>
0 C
,<
O0
DEGRADATION
BY
THE
CROWN-OF-THORNS
J
~
STARFISH
,
oo
I
,
\\
:
:
-",,"
,
•
\
O
I-,
i o
./ t
•
~.~
~11"o
~'"~' ~ o6~',
-J ~ ~ , , ~
I
~
I
,
o
3'~ ~,
-C'
I I
II
I~
I
~=.~
~
~
---
11 t~
o
I
I
~
0
o
~
II , . ~ff~ O'-N
~.-:
ll9
120
R.M.
SEYMOUR
AND
R. H. B R A D B U R Y
the fact that the fine weighting assigns a non-zero weight to the anomalous years 12 and 13. The residuals in both cases were also tested for randomness using the Kendall rank-correlation test (tau) and the turning point test (Kendall, 1973, chapter II), and no significant deviation from randomness was detected at the 5% level. Notice however, from the 95% confidence limits for the data points (Fig. 3) that any such deviation would be more than likely to be swamped by sampling errors. Finally, we note that the values ofp~(0 ~ 1) obtained from these fits do not all lie in the range 0-1. Those outside this range are as follows: for the coarse weighting ( n , p , , ( O ~ 1))=(1, 1.03), (2, 1.08), (3, 1"01); for the fine weighting (n, p ~ ( 0 ~ 1)) = (I, 1.04), (2, 1-09), (3, 1.03). Discussion of these anomalies will be given in section 6.
4. Alternative Models
Though the fit of the main model to the central sector data appears satisfying (Fig. 3), the only real test of its performance is to compare it with possible alternatives. Indeed, it may be thought that the apparent success is simply due to the large number of parameters (nine) involved, and that almost any model with so many parameters would perform just as well, irrespective of biological plausibility. To test this (null) hypothesis, we consider a model of the form (1) with pn(0 ~ 1) = q, as before, but with pn(0 --* 1) a polynomial of degree 4, p,,(O ~ 1) = ao + aim + a2m 2 + a3m 3 + a4m 4
(8)
Taking 1=3, we have nine parameters to fit as before, c = l - q , ao-a4 and ~o, -~t, -~2. The result obtained with the coarse weighting is: c=0.0404, ao=1.515, a1=-2"0859, a2---1"1233, a3=0"0922, a4=-0.1533, .~0=0"7672, .~l = 0-7976, :~2=0.5865 [R2=0-9578; /~2=0.6537]. With the fine weighting: c=0.081, ao=1-1328, a , = - 0 - 4 7 3 4 , a2=-0.6282, a3=0.7664, a4=-0.2304, ~0=0.7213, ~ =0.7922, .~2=0-5126 [R2=0.9512; /~2=0.5802]. The coarse-weighted model is shown in Fig. 4. Clearly, in both cases the fit, as measured by R 2 or/~2, is worse than that of the main model of section 3. However, it is somewhat unclear exactly what this statistic measures (especially when the accuracy of the data is rather uncertain, as in this case). We believe that a better comparison can be made using a likelihood ratio technique (see Appendix). Thus, we are asked to choose between the main model (hypothesis Ho), and the alternative polynomial model (hypothesis Ht) as representations of the true situation, on the basis of the survey evidence over the 24-year period (Table 1). Assuming equal prior probabilities for H0 and H , , the posterior probabilities are as follows: for the coarse-weighted models, P(H0]8)=0.8367, P(H~ [d' ) = 0-1633 (d' = evidence) ; for the fine-weighted models, P ( H o [dr) = 0- 7764, P ( H ~ I g ) =0-2236. Thus, for the coarse weighting, the evidence supports H0 more than five times more strongly than it does H~, and more than three times more strongly in the fine-weighted case. (m=n/lO).
DEGRADATION
BY THE C R O W N - O F - T H O R N S
•
7
r~
I
#
I
B
tq e~
J I
v
STARFISH
i', m
I
A
v
I \
A
v
\
8 0
, Zo/,
f
°
e~
O
O e~ tD
u3
e~ ~:~
.
e~ I I~ I~ O
-
l
to)
°~
Iv i@
i
l i
6
i
i
i
i
o°
~
121
122
R.M.
SEYMOUR
AND
R. H. BRADBURY
In addition, the parameter values obtained with the coarse and fine weightings are not commensurate. Also, the values of p, (0 ~ 1) stray outside the range 0-1 at both ends; for the coarse weighting, the anomalies are: (n, p,(0 ~ 1)) = (0, 1.15), (I, 1.32), (2, 1.14), (21,-0-04), (22,-0.25), (23,-0.51); for the fine weighting they are: (0, 1.13), (1, 1-08), (2, 1.02), (21,-0.01), (22,-0.19), (23,-0.4). However, more biologically meaningful alternative models of the form (1) are possible. In particular, we have already noted in section 3 the anomalies created by the assumption pn(1 ~ 0)=constant; both in the data and the main models this assumption yields values ofpn(0 ~ l) which are too large for small n. From (2) it is clear that these anomalous values would be reduced if c were larger for small n. Thus, as an alternative to the long-term exponential decay in p,(0--* 1), we might suppose such a decay in c = l - p n ( 1 ~ 0 ) . Biologically, this would mean that p,(1 ~ 0) is smaller for small n; i.e. that infestations lasted longer (on average) in the later years than in the earlier years, a feature which would presumably reflect some long-term increase in reef susceptibility (possibly because of increased dominance of reef community structure by fast growing coral species on which Acanthaster feeds preferentially--cf, the discussion in section 6). In view of this, the following forms may be postulated c = 1 - p n ( l -* 0) = c0 exp (-bon)
(9)
with p,(0 --* l) as in (7). This increases the number of parameters from nine to ten. This model was first fitted with b in (7) fixed at 0, thereby reducing the number of parameters back to nine, and forcing all long-term decay into the recovery rate p~ (1 ~ 0). The results obtained are; for the coarse weighting: Co= 0.307, b0 = 0.1782, E=0-724, a=0-2484, k = - 3 . 4 2 1 2 , T= 10.6979, 20=0.6704, 21 =0.8091, 22=0.365 [R2=0"9654; /~2=0"7155; P(HoJS)=0"9321, P(HIJ~')=0.0679]. For the fine weighting: c0=0-303, b0=0 . 1713, E=0-7174, a=0.2619, k = - 3 . 3 4 3 4 , T= 10.7397, 2o=0-6709, 21=0.8469, 22=0.3302 [R2=0.9552; /~2=0-6143; P(H0J4r)=0-9697, P(H~ J~')=0.0303]. The fine-weighted model is shown in Fig. 5. The above models do achieve the aim of bringing p,(0 ~ l) within the range 0-1 (max p,(0 ~ l) = 0.71 and 0-73 for the coarse and fine weightings, respectively, both with n = 3). However, the overall fit is worse than that of the main model as measured by R 2 or/~2, and very considerably worse as measured by P(H0Jd?) and P(HI [g). To see if the exponential decay in c obtained above is genuine, we also fitted the full ten parameter model given by (7) and (9). The results are as follows: for the coarse weightings: c0=0-1289, bo=0-0008, E = 1.2306, b=0-0607, a=0.1976, k = -2.4624, T=11.9998, 2o=0.6772, 2~=0.7907, 22=0.6255 [R2=0.9787; /~2= 0- 7748; P(Ho J°~) = 0.5, P(H~ [ ~') = 0.5]; with the fine weighting: Co= 0.1137, bo = 0.0047, E = 1.2407, b =0-0603, a = 0-2487, k = -2.516, T = 11.9218, .~o= 0.6729, 2, = 0-7919, 22=0.601 [R2=0-9729; /~2=0-7334; P(H0]~')=0-501, P(HjJ~')=0-499]. These results are virtually indistinguishable from those of the nine parameter main model (section 3), with b0 very small in both cases. Thus, when given the choice, the data refuses exponential decay in c but accepts it in p~(0 ~ I). The full ten parameter model also retains the anomalous value ofp,(0 ~ l) found in the main model. It can be concluded that the assumption p~ ( 1 ~ 0) = constant is not sufficient to explain the
DEGRADATION
BY T H E C R O W N - O F - T H O R N S
STARFISH
S
./-.+
,
O
I
I
v
I
#
i
.!-i
7
.Z,/-' : I"
v
t
'
#
+ l
0
t ~
vJ
'
J
:
O
:
l
"1 6
t
l
I
1
-o o
~6~.
123
124
R . M . SEYMOUR AND R. H. BRADBURY
occurrence of these values, and that this hypothesis is reasonable. Further discussion of this matter is given in section 6. 5. The Northern and Southern Sectors
There is insufficient good data available for either the northern or southern sectors of the GBR to permit any sensible fitting of the model (section 3) for these sectors separately. However, we may compare the extant data for these sectors with results from the central sector. Taking the model plot (Fig. 3) as a reasonable overall representation of the central sector data, these comparisons are made in Figs 6 and 7. It seems clear that similar overall processes are at work in these sectors to that in the central sector, but apparently at a lower level of intensity, and with a faster decay towards the end of the period for the northern sector, but slower for the southern. However, there is no striking evidence from these figures of any obvious time-lags which would indicate a process initiated in the northern sector and moving progressively southward through the central and southern sectors in subsequent years. Such a process of build-up in the north followed by a release in a southward moving wave was originally conjectured by Kenchington (1977). To the contrary, the comparisons of Figs 6 and 7 are at least as compatible with the existence of a more or less synchronous process over the entire length of the GBR, but which has maximum strength in the central sector. If there were waves of outbreaks carried from the northern to the central sector by larval transport, then a time-lag of between 2 and 3 years would be expected between corresponding peaks and troughs in the two sectors (larvae produced on outbreaking reefs in the north take this length of time before appearing as adults on central sector reefs). However, the apparent northern sector peak in year 1 does not produce a central sector peak 2 or 3 years later, but is synchronous with such a peak. Similarly for the apparent trough in year 4. Likewise, both processes rise to apparent peaks somewhere between years I4 and 17. Note that the very large northern sector value in year 15 comes from a small sample (nine reefs), and in view of the much more reliable points in years 14 and 16, should be treated with some scepticism. But in any case, the data for the central sector (Table 1) also show a peak for this year. On a more detailed level, the (anomalous) low data point for the central sector in year 19 (Fig. 3) is reflected in the sharp crash from years 18 to 19 in the northern sector. Similarly, the small scale rise and subsequent fall during years 21, 22 and 23 in the central sector, is reflected in a synchronous process in the northern sector. Only year 20 for the central sector seems out of line. Finally, note that the large value in year 13 (four out of five reefs in state l--see Table 2) corresponds to an exactly similar value for the central sector (Table 1). Low confidence should be given to these small sample values; in particular, the value in year 11 is essentially worthless (Table 2). The 0 value in year 9, although good, is too isolated to be easily interpreted. It has to be admitted that all of this is rather impressionistic and open to considerable doubt, due to the inherent uncertainty in the data (section 2). However, if we also take accouint of the southern sector data (Fig. 7), the case for synchronicity
DEGRADATION
BY THE
CROWN-OF-THORNS
STARFISH
125
appears rather stronger. Here, the process looks highly compatible with that of the central sector, though perhaps with greater damping. Certainly, any apparent timelags seem to indicate (if anything) that the southern sector wave precedes that of the central sector. 6. Discussion
It can be concluded from the analyses of sections 3 and 4 that the main model (1) with 1= 3, pn(1 ~ 0 ) = q, a constant, and pn(0 ~ 1) given by (7), yields a satisfactory representation of the (weighted) data for the central sector of the GBR (Fig. 3). Furthermore, commensurate results were obtained with the two weightings considered (coarse and fine--see section 2). Also, this model performed better than some possible rivals (section 4), and in particular, better than the (biologically meaningless) null model in which p~(0 ~ 1) is a polynomial of degree 4 in n. The good fit obtained (and the order of magnitude of the parameters obtained) does not depend strongly on the weighting chosen, provided the weighting chosen gives low values to data points below the 50% confidence level (with ~ = 0-15--see section 2), and in particular, gives a low weighting to the anomalous value in year 13 (Table 1). For example, a weighting which is rather cruder than either of those used here is obtained from a one-sided normal curve with 5% tail at -0-5 (ty= 0.5/1.96), in the manner described in section 2. In spite of giving larger weight to relatively uncertain data points, this weighting produces results commensurate with those obtained in section 3, but only if year 13 is disregarded. Since this point is based on a small sample, which is more than likely to violate the random sampling assumptions used to obtain the 95% confidence limit (Table 1), this circumstance need not trouble us unduly. An additional source of discomfort for the main model was the anomalous values o f p , ( 0 ~ 1) produced in years 1, 2 and 3 (sections 3 and 4). As shown in section 3, these anomalies are inherent in the data if we take pn(1 ~ 0 ) = c o n s t a n t . Further, while the anomalies can be suppressed by forcing c = 1 - p ~ ( l ~ 0) to decay exponentially with n (9), the data firmly reject this option in favour of the (approximate) constancy ofp~(l ~ 0) when given a choice (section 4). We conclude that the structure of the main model is unlikely to be the source of these anomalies, and that they are most likely to arise from sampling errors in the data during these early (and, in terms of outbreaks, very volatile) years. All the results presented here are obtained with latency period 1= 3 years, the reason for this choice being that it is roughly the known time span from the settlement of starfish larvae to the achievement of average adult size. However, the main model was also investigated with l = 1, 2 and 4 years, but only poor results were obtained. In particular, with l = 4 the number of parameters increases from nine to ten, but nevertheless, a less satisfactory fit to the data is obtained. These facts strengthen our conviction that l= 3 captures a significant biological feature. The infestation recovery rate, pn(l ~ 0 ) = q = l - c , given by the main model was of the order, q~0-9 (for either the coarse or fine weighting). This shows that approximately 90% of infested reefs will revert to low starfish densities in any 3-year period,
126
R. M. S E Y M O U R
AND
R. H. B R A D B U R Y
0
.~
_
QO
/
I
t=
~./ ~
:
/ \
I
-~-
t,,.,
-~b
_,,2_
I /
ol,.~
!
o_
.+.j~'Y
H
.++
m
:°
0
t--
4
~ +,,.,.., 0-~ m
<. [
I
(M
,+P:
,
i
II ~+[
~
l
[ I
12) 0
0
<5
0 o
DEGRADATION
BY
THE
CROWN-OF-THORNS
/
#
STARFISH
~
oo
e-
E i
I
m
,*_
I I to
e-, .m
I~o4
'T ,
i 0O
p..
d...... I
E~ o
t.O
/O
I
O
tJ
tJ
t~
t~
t~
r~
o,J
E'or~8 •
U
O 6
O
8
127
128
R. M.
SEYMOUR
AND
R. H , B R A D B U R Y
a figure which is consistent both with field observations (Johnson et al., 1989), and with studies of the average life span of starfish (5-8 years from settlement--Moran, 1986). The major interest however, is in the time-dependent infestation rate, pn(0--* 1)-(7). The period obtained (section 3) is T,~ 12 years. This is particularly satisfying because it is consistent with field data indicating a re-establishment time of approximately 12 years for fast growing coral species (e.g. Acropora sp.) on which Acanthaster is known to feed preferentially (Moran, 1986; Reichelt et al., 1990). Thus, it appears that this food availability cycle is the primary process governing the rate at which reefs become infested. Note that the period between years 11 and 13 (19771979) for which there is no, or only poor data (Table 1), almost certainly does represent a true low in outbreaks (in spite of the anomalous year 13--see above). In fact, the reason little surveying was done during this period is because there were few reports of outbreaks, and the phenomenon was thought to be dying out. We have good reason therefore, to believe that this ~ 12-year cycle is not just artifactual. The exponential decay term in p~(0---, 1) indicates that, on average over the 24year period, progressively fewer of the uninfested reefs were becoming infested. This would seem to show that reefs which have once been infested do not rapidly return to their prior state of susceptibility, even on a time scale long enough for the faster growing coral species to re-establish. It has been suggested by several authors that long-term community degradation may result from starfish damage to long-lived, slow growing species, which are thought to be the main architects of reef community structure (Done et al., 1989; Endean et al., 1989). Our results lend credance to this view. However, the data sequence is insufficiently long to show whether or not this process is reversible (periodic?). The available geological evidence suggests that it might be (Frankel, 1977; Walbran et al., 1988), though the interpretation of this evidence is controversial (Moran, 1986; Reichelt et al., 1990). A mathematical model of such a long-term cycle based on reversible degradation has recently been developed by Seymour (1990). However, it is equally possible that we are witnessing a relatively rapid, and irreversible transition between two stable states (Bradbury et al., 1985)high coral cover, low starfish density to low coral cover, low starfish density, perhaps with an algal dominated community structure--initiated by some unprecedented event, the most likely candidate for which is removal of starfish predators by human activity (Ormond et al., 1990; Reichelt et al., 1990). The main model of section 3 gives a good representation of the data from the central sector (Fig. 3), and indicates the operation of simple processes manifest over the whole region. As far as it is possible to judge from the meagre data available, similar processes appear to operate in the northern and southern sectors also (section 5). However, there is little real indication of significant time-lags between the three sectors which would substantiate the northern sector outbreak hypothesis of Kenchington (1977)--although possible sampling errors prevent any definitive assertion. Instead, global pulses over the entire GBR are suggested, with maximum strength in the central sector. We believe this is the first (empirical) indication of the existence of such synchronous global pulses of starfish outbreaks. The only other work of which we are aware which explores even the possibility of this phenomenon,
DEGRADATION
BY THE CROWN-OF-THORNS
STARFISH
129
is the series of studies of simple cellular automata models of the GBR by Bradbury et al. (1990) and van der Laan & Bradbury (1990). These authors were able to generate a range of global outbreak behaviours in their models by varying connectivity between reefs and the intensity of fish predation on starfish. In particular, they generated patterns of many small, isolated, non-propagating outbreaks at low levels of connectivity or at high levels of fish predation. This is the pattern which some authors (e.g. Cameron, 1977) believe to be normal for the system. Travelling waves of outbreaks were generated broadly under intermediate levels of connectivity or predation, while global pulses were generated under conditions of high connectivity and/or low predation. Pulses seem to be the upper bound to the global dynamics of the system, and so suggest that some process is operating close to or at a critical limit (cf. Seymour, 1989, 1990). Since connectivity is presumably a permanent feature of the system, it is perhaps more likely that this is predation. Finally, the main model of section 3 can be used to predict the behaviour of outbreaks over the central sector of the GBR. The prediction to the year 2000 is
1.0
0.5
0
I
2
3
0'11
0"15
I n 12 x
0"08
4
5
I 6
I 7
8
9
I0
I
I415 6 I, 819 °l I 0"2
0'23
0"23
0'21
0"18
0'14
0'1
0"07'
FIG. 8. Prediction of proportion of reefs in state ! (x) up to year 2000; n = y e a r ( 0 = 1989). All points are 4-0.13 (empirical 95% confidence limits--see text).
130
R. M. S E Y M O U R A N D R. H. B R A D B U R Y TABLE 5 The proportion x = s/Y~ o f reefs in state 1 in recent surveys
Year
1990/91
1991/92
Northern sector Central sector Southern sector
0.0 (E=34) 0-16 (Z = 43) 0-02 (Z = 44)
- - (Z=0) 0.10 (Z = 31 ) 0.33 (E = 3)
E is the number of reefs surveyed in each sector.
given in Fig. 8. N o t e t h a t o u t b r e a k s a g a i n reach a p e a k in the mid-1990s. T h e 95% c o n f i d e n c e limits are o b t a i n e d using the e m p i r i c a l s t a n d a r d d e v i a t i o n o b t a i n e d f r o m the d a t a o f T a b l e I ( K e n d a l l , 1973, c h a p t e r 9). O f course, this p r e d i c t i o n is p r e d i c a t e d on the a s s u m p t i o n t h a t the l o n g - t e r m d e g r a d a t i o n o f reef c o m m u n i t y structure, implicit in the e x p o n e n t i a l d e c a y t e r m in the m o d e l , will continue. T o verify o r refute this, we suggest that survey d a t a on starfish a b u n d a n c e c o n t i n u e to be collected, a n d t h a t d e t a i l e d fieldwork o n the effects o f starfish o n long-lived coral species s h o u l d be c a r r i e d out, with a view to i l l u m i n a t i n g the r e - e s t a b l i s h m e n t c a p a c i t y (if a n y ) o f reef c o m m u n i t i e s o v e r the l o n g e r term. W h e t h e r o r n o t such r e c o v e r y is p o s s i b l e is, in o u r view, the key to a n y r a t i o n a l m a n a g e m e n t p o l i c y ( B r a d b u r y & Reichelt, 1982). Note added in proof : Since the manuscript was completed, Dr Peter Moran of the Australian Institute of Marine Science has kindly made available a summary of the unpublished results of his latest surveys of starfish outbreaks. The organization of these surveys has been redesigned recently, so that each designated survey area of the GBR is sampled once during each ist July to 30th June year. This means that the new sampling year is approximately 6 months out of phase with the (calendar) year on which previous data were based. The data at hand concern the finalized surveys for the 1990/91 year, and the as yet incomplete surveys for the 1991/92 year. Fortunately the central sector has been well sampled in this year. The data are summarized in Table 5. By plotting the appropriate data points at 1990.5 (= 1.5) and 1991.5 (=2.5) on Fig. 8, it can be seen that the model prediction here is accurate to within the 95% confidence limit for these 2 years.
REFERENCES BECKER, N. G. (1989). Analysis of Infectious Disease Data. Monographson Statistics and Applied Probability. pp. 221. London: Chapman and Hall. BroKEr.AND, C. (1989). The Faustian traits of the crown-of-thorns starfish. Am. Sci. 77, 154-163. BRADBURV, R. H., HAMMOND,L. S., MORAN, P. J. & REICHEt.'r, R. E. (1985). Coral reef communities and the crown-of-thorns starfish: evidence for qualitatively stable cycles. J. theor. Biol. 113, 69-80. BRADBURV, R. H. & MUNDY, C. N. (1989). Biomass Yields and Geography of Large Marine Ecosystems. (Sherman, K. & Alexander, L. M., eds) pp. 43-167. Boulder, CO: Westview. BRAOBURV, R. H. & REICHELT, R. E. (1982). The reef and man: rationalizing management through ecological theory. Proc. IVth Int. Coral Reef Syrup. 1,219-223. BRADBURV, R. H., VAN t)ER LAArq,J. D. & MAcDONALO, D. (1990). Modelling the effects of predation and dispersal on the generation of waves of starfish outbreaks. Math. Comput. Model. 13, 61-67. CAMERON, A. M. (1977). Acanthaster and coral reefs: population outbreaks of a rare and specialised carnivore in a complex high-diversity system. Proc. Hh'd. Int. Coral ReefSymp. 1, 193-199. DONE, T. J., OSBORNE, K. & NAVtN, K. F. (1989). Recovery of corals post-Acanthaster: progress and prospects. Prac. Vlth Int. Coral Reef Syrup. 2, 137-142.
DEGRADATION
BY T H E C R O W N - O F - T H O R N S
STARFISH
131
ENDEAN, R., CAMERON, A. M. 4~: DEVANTIER, L. M. (1989). Acanthasterplanci predation on massive
corals: the myth of rapid recovery of devastated reefs. Proc. Vlth Int. Coral Reef Syrup. 2, 143-148. FRANKEL, E. (1977). Previous Acanthaster aggregations in the Great Barrier Reef. Proc. lllrd Int. Coral
Reef Syrup. 1,201-208. FRAUENTHAL,J. G. (1980). Mathematical Modelling in Epidemiology. pp. 118. Berlin: Springer Verlag. HOWSON, C. & URBACH, P. (1991). Bayesian reasoning in science. Nature, Lond. 350, 371-37:4. JOHNSON, D. B., BAss, D. K., MILLER-SMITH, B. A., MORAN, P. J., MUNDY, C. N. & SPEARE, P. J. (1989). Outbreaks of the crown-of-thorns starfish (Acanthaster planci) on the Great Barrier Reef: Results of surveys 1986-1988. Proc. VIth Int. Coral Reef Syrup. 2, 165-169. KENCHINGTON, R. A. K. (1977). Growth and recruitment ofAcanthasterplanci (L.) on the Great Barrier Reef. Biol. Conserv. 11, 103-118. KENDALL, M. G. (1973). Time Series. pp. 197. London: Griffin. LAAN, J. D. VAN DER & BRADaURV, R. H. (1990). Futures for the Great Barrier Reef ecosystem. Math. Comput. Model. 14, 705-709. MORAN, P. J. (1986). The Acanthaster phenomenon. Oceanogr. Mar. Biol. Ann. ReD. 24, 379--480. MORAN, P. J., REiCHELT, R. E. & BRADnURY, R. H. (1986). Assessment of the geological evidence for previous Acanthaster outbreaks. Coral Reefs 4, 235-238. ORMOND, R., BRADBURY, R. H., BAINBRIDGE, S., FABRICIUS, K., KEESlNG, J., DEVANTIER, L., MEDLV, P. & STEVEN, A. (1990). Test of a model of regulation of crown-of-thorns starfish by fish predators. In: The Acanthaster Phenomenon: A Modelling Approach. Lecture Notes in Biomath., no. 88 (Bradbury, R. H., ed.) pp. 189-207. Heidelberg: Springer Verlag. REICHELT, R. E., BRDBURY, R. H. & MORAN, P. J. (1990). The crown-of-thorns starfish, Acanthaster planci, on the Great Barrier Reef. Math. Comput. Model. 13(6), 45-60. REICUELT, R. E., GREVE, W., BRADatmY, R. H. & MORAN, P. J. (1990). Acanthaster planci on the Great Barrier Reef: a starfish-coral site model. Ecol. Model. 49, 153-177. SEYMOOR, R. M. (1989). In Acanthasterplanci a near-optimal predator? Ecol. Model. 46, 239-260. SEYMOUR, R. M. (1990). Very long period cycles in a near optimal model of the population dynamics of Acanthaster planci. IMA J. Math. AppM. Med. Biol. 7, 157-174. WALBRAN, P. D., HENDERSON, R. A., JULL, A. J. T. & HEAD, M. J. (1988). Evidence from sediments of long-term Acanthaster planci predation on corals of the Great Barrier Reef. Science 245, 847-850.
APPENDIX
Let N be the (very large) number of reefs in a given sector (N is in the region 8001000), and let a be the proportion of these which are infected (state 1) in a given year. Assuming that the infected reefs are distributed at random, then a random sample of n reefs (n<
I"
a0--
•
(A.D
N o w , f r o m the result o f a given s u r v e y (i.e. given r a n d n), we wish to o b t a i n i n f o r m a t i o n a b o u t a. L e t He be the h y p o t h e s i s t h a t a E [ ~ , ~ + d ~ ) ; ( 0 < ~ < 1). T h e e v i d e n c e f o r this h y p o t h e s i s is the s u r v e y result, d r = {r, n}. I t follows f r o m ( A . I ) t h a t the c o n d i t i o n a l p r o b a b i l i t y , P ( 8 [ H ¢ ) , is P(#,H~)=(:)~r(1-
~)n-r+ ~(d~).
(A.2)
H o w e v e r , it is the p r o b a b i l i t y o f H~ c o n d i t i o n a l o n the e v i d e n c e ; i . e . P ( H g [ $ ' ) t h a t is o f a c t u a l interest. T h i s c a n be o b t a i n e d f r o m B a y e s ' s t h e o r e m P(H~I¢) = P(¢IHDP(H~)/P(8).
Therefore only the prior probabilities
P(H~) and
P ( 8 ) need to be assigned.
(A.3)
132
R.M.
SEYMOUR
AND
R. H. B R A D B U R Y
In the absence of any evidence, the simplest and most natural choice is to assume that every value of a is equally likely; i.e. that the distribution is uniform,
P(H~) = d4.
(A.4)
Similarly, for a given survey of n reefs, there are n + 1 possible values of r (0 to n), and prior to the result of an actual survey, it is most natural to assume that each of these is equally likely; i.e.
e(8) = I/(n+ l);
dr= {r, n}.
(A.5)
It now follows from (A.2)-(A.5) that the distribution we are seeking is given by
[Note that this is indeed a probability distribution; i.e. ~J0P(tI~I~)= 11. It is now possible to obtain confidence estimates using (A.6)--and therefore based on the random sampling and uniformity assumptions described above, and which define a (null) prior hypothesis. Thus, given ~' = {r, n}, and an error level 4, f rln+ ~
e(r/n - 4 < a < r/n + ~)
P(H~ Ig).
(A.7)
= ~ r/n--
To find 95% confidence limits, it is required to find 4o and it with 0 < ~ o < ~ j < 1, satisfying
f¢
" P(H, Ig ) =0.95.
(A.8)
o
If r = 0 , then 40=0 and 4~ = 1 - a ~ ( a =0.05, u = 1 / ( n + 1)), while if r = n , 4 o = a ~ and 4w= 1. On the other hand, i f 0 < r < n , we choose those values of 40 and ~j for which the length of the confidence interval is a minimum. The condition for this is easily seen to be ~(1 - ~ o ) # - ' = ~ ( 1
- 4 , ) n-r.
(A.9)
Since (A.6) has a unique maximum at 4 = r/n, it follows from (A.9) that 40 < r/n < 41. Given r and n, (A.8) and (A.9) determine 40 and 4t uniquely. Finally, the likelihood ratio procedure for distinguishing between two model hypotheses used in section 4 is considered. Given two models M0, M~ of the form (1)--with parameters specified--we wish to determine which of them best represents the data. Since the parameters have already been determined by least squares estimation from the data, M0 and M~ can be regarded as the best adapted models of their type, so that the remaining choice is between types. These models generate sequences of points, M,- = {~n [0 <_n < 23 } (i = 0, 1), where ~in is the proportion of reefs in state 1 in year n according to M,.. It follows that, if 8 = {(sn, Z , ) 1 0 < n < 2 3 } i s the survey
DEGRADATION
BY T H E C R O W N - O F - T H O R N S
STARFISH
133
evidence (Table 1), and H~ is the hypothesis that Mt is a representation of the true situation, then [see (A.I)]
p(SlH,) = ii(Zs).~,.( 1 _ ~,)~-s
(A.10)
where the product is taken over all years n for which data exist. Define the likelihood ratio (of Hi over H0) to be L R = P ( ~ I H t ) / P ( 8 t H o ). It now follows from Bayes's theorem that
e(Ho [~) = P(Ho)/(P(H,)LR + P(Ho))
(A. 1l)
and P(Ht [g) = 1 - P(Ho[$'), where P(Ho), P(H)) are prior probabilities. Clearly, our prior assumption must be that H0 and/-/i are equally likely, so that P(Ho)= P ( H 0 = 1/2. It then follows from (A.11) that
e(Ho]8) = I /(LR + l).
(A.12)
It is (A.10) and (A. 12) which are used to compute the posterior probabilities given in section 4.