THE U N I V E R S I T Y O F M I C IS I G A N COLLEGE OF LITERATURE, SCIENCE, AND THE ARTS Computer and Communication Sciences Department MULTI-MODEL MULTI-FORMALISM MODE LLI NG: AN ECOSYSTEM EXAMPLE Bernard P. Zeigler August 1977 Logic of Computers Group Computer and Communication Sciences Department Technical Report No. 206 with assistance from: National Science Foundation Grant No. MCS76-04297

Multi-model Multi-f ormalism Modelling: An Ecosystem Example by Bernard P. Zeigler Department of Applied Mathematics The Weizmann Institute of Science Rehovot, Israel and Logic of Computers Group Department of Computer and Communication Sciences The University of Michigan Ann Arbor, Michigan Prepared for "Theoretical Ecosystems", Ed. E. Halfon, Academic Press (In preparation) This research was supported in part by National Science Foundation Grant No. MCS76-04297

1. 1. Introduction Ecosystems, as examples of large scale multifaceted systems, require that: a multiplicity of models be developed since a single all encompassiug model., however desirable as a conceptual goal, is not a practical object. By decomposing questions and modelling objectives into an ordered structure of elements called experimental frames (Zeigler, 1976a), useful partial models may be constructed, validated and employed, each one attuned to a partictilar experimental frame. Concormmitant with the pluralism of such partial modecls is the recognition that models are expressible in different formalisms, each offering conceptual and computational advantages within its domain of application (Zeigler and Barto, 1977). But now, in addition to the familiar activities involving construction and validation of individual models, there is required a host of organizational activities aimed at integrating the collection of models into a synergistic whole. Our belief is that the computer can be programmed to aid in executing these activities to a much greater degree than it is doing today. We have recently sketched a theoretical basis for structuring the organization of partial models (Zeigler, 1977a). In this paper, we illustrate our approach by considering an ecosystem example in some detail. After briefly reviewing this approach, we discuss its application to the patch structured universes employed by Huffaker (1958, 1963) to study predator-prey coexistence. We show how the approach facilitates the development of mutually supportive detailed and abstract models that in conjunction, provide both accurate ecological realism at one extreme and general insight into the essential mechanisms at work, at the other.

2. The Ecosystem: Questions of Interest and Models The problem of predator-prey coexistence in patchy environments has received muc:h theoretical attention of late (Levin 1976, lass;ell and May, i. //I, Maynard Smith, 1974). Most of this work has employed the conventional differential and difference equation formalism but following on suggestiJonx.s stressing the importance of discrete processes (Maynard Smith, 1974), we have shown (Zeigler, 1977b) that the discrete event formalism and associated simulation languages can provide effective comprehensible exp1laation.P predator-prey co-existence. Also. there have been very few attempts to fit the theoretical models to labor-oary or field data. In contrast, our general approach is illustrated in this case by our development of four related models, expressed in either the differential equation or discrete event formalism, and constructed at different levels of abstraction, ranging from the most detailed level where close comparison of model behavior with experimental data is possible, to the most abstract where overall properties are discernable in a relatively simple manner. The real system to which the modelling is directly addressed is that of the controlled universes constructed by Huffaker (1958, 1962). Basically, tiese consist of spatial arrays of oranges (patches) of controllable nutritional valu, inhabitable by prey and predator mitesand interconnected by migration pathways of controllable difficulty. Many other discrete food unit environments fit this general form. Our most detailed model, the base model, is of the stochastic differential equation type. In it, the local state (situation on each patch) is determined by a Lotka-Volterra type differential equation governing the joint food, prey

and predator dynamics; the impetus for emigration and the effect of immigration are logically determined from the local state (food, prey, predator); and the migration process is of the stochastic random walk 'a-ricl y. This model enables us to identify the parameters of the local Lotka-Volterra dynamics from data for single patches. It is not feasible for computer simulation however, and this motivates the construction of a second model which is both simulateable and amenable to validation against data collected from universes in which the effect of migration is at issue. This (second) model is of the stochastic discrete event type and is simulated in SIMSCRIPT, a well known discrete event language. The model keeps track of the same state variables as its predecessor but updates them only at "event times". The tables required for scheduling events and executing the updates were derived by appropriately partitioning the local state space and summarizing the Lotka-Volterra trajectories between partition boundaries. (This required a once-and-for-all simulation of the Lotka-Volterra equations.) This technique for representing differential equation models in summary form as discrete event models is quite general (Zeigler, 1977c). Our third model in the hierarchy is also of the discrete event type. It is an astraction of its predecessor, in which the local situation is represented by a small number of discrete states [empty, some prey, etc.] and the migration processes are also suitably simplified. Since actual population counts have been discarded, this model cannot make quantitative global population predictions. The model lends itself however to convenient parameter study of persistence and the development of patterned interaction.

4. Our last lumped model is of the deterministic different equation type. It describes the global behavior of its predecessor operating under the so called "random phase-random space" mode (Zeigler, 197Lb). The model yields simple algebraic expressions for the equilibrium distribution of patch states and thus explains the form of the dependence of persistence oiL migration parameters. 3. Organization of Questions and Models The integration and organiztation of the above models is achieved within the formal system sugges.ed by Zeigler (1976a, 1977b). The following is an informal review of the concepts involved. We distinguish the following elements. -- a collection of experimental frames. A frame E C represents a restricted set of questions by specifying the restrictions on experimental access to the real system sufficient to answer them. Such a frame E dete-mtnies a collection of data sets ~ (E), such that each;) E (E) is an a priori possible result of complete data acquisition within -rme E 4a -- the real system, is comprised of the specific data that has been, or would be, collected by experimenting with the system. Thus R associates with each experimental frame E a unique data set A (E) EC ((E) &/^ -- The doinain of possible models. These are assumed to be transition systems which are specifiable at various levels of structure and behavior and within various short-hand conventions such as the sequential machines, discrete event and differential equation formalisms.

A full description of these basic eleuments and the concepts they embody may be found in Zeigler, 1976, Chapters 2 and 11. In a moment, we shall formulate these elements in the context of patch structured universes and the experiments of Huffaker (1958, 1962) in particular. Roughly, the "experimental frames" will encode the various choices of observables (species counts in patches) and conditions (initial stocking of species, sLtucture of universe) undcer which experiments were run. The "real system" is the data collectarblc by making the implied observations under the given conditions. Finally the "models" are the various distributed and lumped models which can be postulated to account for the obse-rved data and to predict thee results of future experiments.

We can imagine an ideal state of affairs in which for each frame E C there is a known model H C6< which "best" answers the questions posable in E. By "best" we mean that the model can reproduce without erro'r the data set ~(E) in a manner which requires the least consumption of computer resources. Realistically, this ideal is not realizable after a necessarily finite span of data acquisition. The dynamics of modelling concern successive approximations to the ideal. We formulate the problem as follows: At any time t, the dat:a already acquired in frame E will be some subset (E)g (E). Many, perhaps most frames, will not even havze )een considered. Those that have, form tlJe subset t = {EI t (E) 7~. Similarly, only a small subset. t of the possible models _ will have been constructed as potential model candidates. Thus the state of affairs at any time t is reflected in the triple 4. Experimental Framest Let us examine the elements of the triple at time t = June 1961, the date of the last observation recorded by Huffaker (1962). There are four main types of experimental frames. As displayed in Table 1, these types are distinguished by the descriptors "global", "local", "total", and "occupancy". The "global" descriptor refers to the fact that all cells (locations where an orange or a substitute rubber ball may be placed) in the universe are being observed. In contrast, in the "local" condition, only some subset of the cells are of interest. The "total" descriptor refers to the fact the quantities of interest in a frame have been totalled to produce aggregate quantities, so that only those aggregates are observable in the frame. acinally, the "occupancy" descriptor refers to the fact that a

frame permits only the observation of discrete occuplancy states, such as whether or not a cell is empty, whether or not a prey colony has been established on the cell, and whether there are no, few, or many predators present. Table 1 also summarizes the kinds of questiorns associated with each frame. The "occupancy" frames are the most restricted, Nonetheless they permit consideration of persistence of predator —prey relations since to determine whether there are any prey or predators requires only a binary categorization (present/not present) for each cello. At the other extre>.me, the "global" frames permit observation of detailed spatial distribution off species. The "total" frames correspond to classical populations in which spatial structure has been averaged out. It is evident that certain frames are potentially more informative than others. In a moment, we shall formally characterize this fact in terms of the "derivability" relation (Zeigler, 1977a). The experimental frames 8 are defined in Table 2. Each frame names a set of variables of interest, called the compare variables, and a set of variables determining the conditions under which experiments are to be performed, called the control variables. The most inclusive frame, global Efoodpreypred specifies as compare variables: food amount, prey, and food Prey,pred predator population counts in each cell. This constitutes a total of 3N variables where N is the number of cells in the universe. There are no control variables for this frame. An example

aDao. Ex-~erimental 'rames and Their Posable uestions Plane of Frame Jescriptiotrq Posable Questions Concern: (representative frames) given Eglobala food amount and population Spatial characteristics of variables for each cell predator-prey, prey-food interaction food amount and population variables for a subset of cells Predator-prey, prey-food, (degree of abstraction O) keeping all others zero interaction in local patch globalitztal fgood amount and population Space averaged population sizes in variables totalled over all cells predator-prey, balance of prey-food, (degree of abstraction I) interaction (classical lumped (degree of abstraction = 1) populations) Ecccupancy,global Discrete food and population Persistence of predator-prey, states for each cell (empty, some balance of prey-food, interactions; (degree of abstraction I= I prey, many prey, some predators, effect of cell geometry etc.) Eoccupancy,%lobal,total Totals of cells in the various Persistence of predator-prey, balance of prey-food, interactions states as given in under random phase-random (degree of abstraction = 2) Eoccupancyglobal space -or(ditlans.

T a.,_ 2 xper~er.zl Fraes and Releva- Dat~ao Sets in Huffaker Universe Frame, -E Jescristion Associated Data SetulP (E) (Table 3) Eglobai Food, prey, predator variables missing food,prey,pred for each cell global Efoo doprey/prd MFood, prey variables for each fdprey/pred cell in absence of predator missing.global Food variable for each cell missing (but partial descriptions of rood/prey,pred in absence of prey and predator the orange replenishment schedules used are given) Eglobal,total Totals of food, prey and 58, II(A-I), Figs. 9-18 Efood,prey,pred predator over all cells 63, 11(3,4), Figs. 3,4 63 I-(4,5,6), Fig. 5 Eglobal,total Totals of food, prey over all 58, I(A,B,C) Figs. 6, 7, 8 food,prey/pred cells in absence of predator 63, E-2, Fig. 2 Eglobal,total Total of food over all cells food/prey,pred in the absence of prey and missing predator -.-.- -S..

Table 2 (continuad) Frame Description Associated Data Sets HoQ(E) (Table 3) local foodpreypred Focd,pray,predator variables for a subset of cells keeping missing all others zero Elocal food,prey/pred Food, prey variables for a missing subset of cells in absence of predator Local fodd/prey,pred Food variable for a subset of missing (some very incomplete cells in absence of prey and description of orange quality predator and spoilage rates given) local,total Elocalodptota d Totals of food, prey, predator 58, II A, Fig. 9 variables for a subset of cells 58, II B, Fig. 10 keeping all others zero 58, II C, Fig. 11 local, total food;prey/pred Totals of food, prey variables for a 58, IA, Fig. 6; Initial Parts of: subset of cells in absence of 58, TB, Fig. 7 predator, keeping all others zero 58, IC, Fig. 8 63, E2, Fig. 2 Eoccupancy, global E ccupancyglobal Occupancy states of prey and predAtor 58, II I, iig. 18 for each cell 63, TT-(3,4), Figs. 3,4 Eoccupancy, global, total. E ccupancyglobaltotal Totals of cells in prey and predator computable from: occupancy states 58, II I, Fig. 18 63, II-(3,4), Figs. 3,4

To ble 3 Data Elements For Huffaker Universes ~(t 1961) Data Element Key* Description (as.. _iven, by luffaker) 58, IA, Fig. 6 Predators Absent, Sltplest Universe, Four Iarge Anreas of Food, Grouped at Adjacent, Joined Positions 58, IB, Fig. 7 Predators Absent, Four Large Areas of Food Widely Dispersed 58, IC, Fig. 8 Predators Absent, 20 Small Areas of Food Alternasting With 20 Positions With No Food 58, IIA, Fig. 9 Predators Present, Simplest Universe, Four Large Areas of Food, Grouped at Adjacent Joined Positions 589 IIB, Fig. 10 Predators Present, Eight Large Areas of Food, GCxwuped at Adjacent Joined Positions 58, IIC, Fig. 11 Predators Present, Six I:hole Oranges as Food, Grouped at Adjacent Joined Positions 58, IID, Fig. 12 Predators Present, Four Large Areasof Food Widely Dispersei 58, IIE, Fig. 13 Predators Present, Eight Large Areas of Food Widely Dispersed 58, IIF, Figs. 14, 15 Predators Present, 20 Small Areas of Food Alternating with 20 Foodless Positions 58, IIG, Fig. 16 Predators Present, 40 Small Areas of Food Occupyilo-ng All Positions 58, IIH, Fig. 17 Predators Present, 120 Small Areas of Food Occypying Al]. 120 Positions (Barriers to Migration Added) 58, II I, Fig. 18 Predator-Prey Oscillations, 120 Small Areas of Food Occupying all 120 Positions (Barriers to Migration Added) 63, E-2, Fig. 2 Predators Absent, Complex 3-Shelf Universe, 210 Small Areas of Food 63,II-(3,4), Figs. 3,4 Predators Present, 3-Shelf Universe, 252 Small Areas of Fo 63, I-(4,5,6), Fig. 5 Predators Present, Complex 3-Shelf Universe, 252 Larger Areas of Food * (58, IA, Fig. 6) denotes that the data set is presented in Fig. 6 and discussed in Section IA of Huffaker (1958).

8. of a frame which has a non empty set of control variables is Ef o food,prey/p red whose compare variables are all 2N food and prey variables and whose control variables are all N predator variables. In a frame with no control variables, compare variables readings are recorded against time for the duration of any particular experiment (See Fig. 1). This yields a time function (aloo called a segment, or trajectory) which we refer to as a data element belonging to the frame. The set of all. such data elements observable in a particular Hluffaket Universe is the dalta set (E) assigned by such a real system to frame E. The set of all such possible data sets assignab by the possible Huffaker Universe is i (;1). In Table 3, we have listed the data elements recorded by Huffaker (1958, 1962). The collection of these data elements constitutes the real system data till time t = June 1961. In Table 2, we distribute these data elements among the experimental frames. The data elements assciated with a frame E in Table 2 constitute the subset (E) ofx (E), nalnel.y, thcdata acquired until time t in frame E. Table 2 displays some frames as having "missing" data sets. A frame such as Efglobal for which this is true is conceivable, i.e., it is food,prey,pred an element of, but up to time t = 1961, no data has been collected for il i.e., it is not an element of &(t = 1961) global food,prey,pred "missing" because the food amount at each individual orange is not recorded in the Huffaker experiments, even though certain aggregate utiliza are. Although many of the conceivable frames were actually realized in the Huffaker experiments, in the current modelling effort we found that the missing frames often contained information which could have been extremely helpful. tne of the benefits of representation of experimentation in the

experimental frame fonrmalism is that certain experiments may be suggested by the logical structure of frame organization which may turn out to be crucial in later modelling. These mnight not have been thought o}f-: i n an-x unstructured experimental approach. In a frame with control variables specified, the above concepts hold except that the compare variable readings are recorded only so long as the control variables remain zero. Thus for example in the class of frames local denoted by Eoca d, each frame specifies a subset of cells S sutch food,prey pred' that the food amount, prey, and predator counts of cells in S are the compare variables, and all other prey and predator densities are control variables. Data is collectable within such a frame so long as no prey or predators establish themselves on oranges other than those in S. When S consists of a single orange, such data give a picture of the local interaction of food, prey and predators uncontaminated by colony establishment on other oranges,

,/U U 0 0 / a - U U 0/o o00 0oo o o /o o o o so oo o o / a) (Captions to!b' found atL cnd o: p-aper) C 5000 0.>_ 4000, 0 o. 3000 J 2000 o M b) ' ~000 E z Figure 1

10. or subsequent remigration from these colonies. Trh:i[le Huf faker did no experiments with single oranges, the saLie principle holds when S is taken as the subset of initially seeded oranges, given ITu~f~ &!,r's observation that migration occurs only due to food depletion or ovec — population (indeed, this is the basis of our discrete event models). Fo) — example, see Fig. 1 We note that "keeping control variables zero" is a special case of tlhv "range of validity" specification given by Zeigler (1976, Chap. II). 4.1. Organization of Frames The frames in 6 are partially ordered by a relation "is derivable from" or in short" t ". E' g E means that the restrictions on data acquisition imposed in frame E' are over and above those in frame E. As a consequence, data collectable in E' can also be deduced from data. collectable in E and questions posable in E' are posable in E as well. Formalizing one step further, we require an onto mapping from (E ) to 3(E') where mapping D to D' has the interpretation that data set D' is derivable from data set D (we write D' < D) by employing the unique set of operating specified by the pair (E', E). Such operations will in general be information destroying in nature, so that, questions answerable given D cannot be answered given D' Figure 2 displays the "is derivable from" relation in our current example. Three types of operations are employed in this case to reduce data sets, one to another. This is apparent in the following definition given first for frames having no control variables:

gl. obal food, prey, pred _local f o)d p pedy, prEld.d/edp -rey,pre d Eglobal, total Eoccupancy,global food/preypred food,prey,pred occupancy Efood/prey,pred occupancy,global, total rfood,prey,pred Figure 2

11. E < E' if: 1. the compare variables of E are a subset of compare variables of E' or 2. the compare variables of E are simple sums of the compare variables of E', or 3. the compare variables of E are obtainable by discretizing the compare variables of E' (for example, the variable with range {empty,some prey, maximum prey} is obtainable from the variable "prey count", or 4. any composition of the above The operation types associated with 1., 2., 3. above are selection, aggregationand coarsening, respectively.

12. For frames having control variables thle definition is E E' ifa) the compare variables of E may be computed from the comparevariables of E' employing selection s aggregation or coarsening, in aDny composition (as above). b) the control variables of F include those of E', and if there are any additional control vari,!)les in E, they may be comrputed fromn the compare variables of E' employing selection, aggregation and coarsen. — ing, in any composition. In Figure 2 and Table 2, the frames are further organized into pranes. Frames in the same plane are relatable using only the selection operation. If E < E', then we place E on a lower plane than E' if at least one aggregation or coarsening operation must be used to derive the variables of E from those of E'. In fact, the minimal number of such operations necessary to make this derivation is a measure of the distance between planes. In particular, the distance from the base, or most inclusive, plane to a given plane is a measure of the degree of abstraction embodied by the latter plane.

13. Note however that the < relation is a partial order so that there may be more than one distinct plane with the same degree of abstraction. Table 1 shows the naturalness of tbie plane notion in the ecosystem context by relating the plane of a frame to the kind of question posabl:.c with in it. As expected, the higher the degree of abstraction, th.!e mlo)rnC restricted the questions posable. But note that the global/local distilnction does not involve a change in degree of abstraction. 5. Constructed Models Table 4 provides summary descriptions of the models constructed to date, ) t ].977, The conceptual basis for the Occupancy and Random-Phase-Space models has been fully described by Zeigler (1977b). We proceed to describe our base model and the discrete event lumped model derived from it. 5.1. Base Model The base model postulates food, prey, predator interaction on an orange in isolation to be specified by the differential equation: dr _ t -ux pos(r)....la) dx= (b pos(r) - d)x - cxy....lb) dt dt where r is the food amount (measured in fraction of unused orange surface), x, y are prey and predator population sizes, and pos(r) is 1 if r is positive and O otherwise, The meaning attached to the 6 parameters involved is given Table 5. The underlying time unit is one day.

Table 4 Models and Brief DescrLption Model Description (Formalism) Base (Combined differential Local state (situation on each orange) detelrmined equation-discrete by Lotka-Volterra type differential equation event; Combined stochastic- governing joint food, prey, and predator dynamlics' deterministic) impetus for emigration an-d the effect: of immigration are logically determined from local. state (food, prey, predator); migration process J. of stochastic random walk variety; orange replacement schedule simulates that emipl.oyed by Buiffaker. Lumpre Keeps track of same state variables as base rnodel (Discrete Event; Combined ": updates them only at "event." times. The t abl) Stochastic-detenninistic) required fori scheduling events and executing the updates were derived from base model local interaction as discussed in text. Occupancy The local. situation is represented by a smal] (Discrete Event; Combined number of discrete states (empty, some prey, etcStochastic-Deterministic) and the migration processes are simplified to the Bernoul-i trial type. Random Phase-Space (Differential Equation; Derived from Occupancy Model under random phlase Deterministic) -random space hyothesis. Describes the dynamics the occupancy probabilities of the discrete st:ate, when occupancy model is operating in rafidom phase-random space mode.

Table 5 Use of Experimental Frames in Identifyin lModel Parameters Model Component Parameter Description Identified in Experimental V)r-?e:7.e Local Interaction local,total (large population b prey birth ratc food, prey/pred model) d prey death rate u prey food uitilization predator death rate local,total foodpreyprey c predation rate H c t predation eficiency 0 Food Replenishment threshold threshold on prey population Elocaltotal below which orange ijs replaced Prey migration pyrem prey fraction remaining Eglobal total after emigration food,prey/predt pysurvive probability of migrating prey finding a cell meanpysearch mean search time for prey finding a cell pydifl,pydif2 prey diffusivities in horizontal and vertical direction Pred Migration globaltotal pdrem analogous Efood prey proCe pdsurvive to meanpdsearch prey pddifl,pddif2 Local Interaction d" predator death probability Eglobal,total (small population at low prey size food,prey,pred model) c" minimum prey required to global,total initiate predator reproduction food,prey,pred c fraction of prey used to create _global,total 1 predator Efood,prey,pred

14. Equation la asserts that food utilization is proportional to prey density (recipient limited i nteractior) so long as food remain,.;. Eqs. la and lb are Lotka-Volte(rrra rel.auiIoriri withollt;el. —limitat:ior. We assume, following HIuffaker's observations, that prey migrate only Wt,'eliw food is exhausted (r first becomes 0). The migration is effected as follonrs;: 1. When (and if ) food is exhausted, a fraction pyrem of the current prey rel,?ain on the orange (and are subject to the dynamics of Eq. 1). 2) Of the migrating prey [(l-pyr.._m) tiices current population)], a fractions pysurvive are assumed to actualiy reach a cell (the rest are lost to the system) 3. For each of the migrating individuals, a search time T is sampled from an exponential distribution with mean imeanpzsearch. 4. The cell assigned to the individual is computed by quantizing spatial coordinates derived from normal distributions (indlependent for each dimension) with mean, the current cell location and standard deviation pydit (in case of horizontal dimensions) or pydif2 2T (in case of the vertical dimension). 5. After time T has elapsed, the individual is added to the population of his assigned cell if food remains there; otherwise, with probability pysurvive he is sent to step 3) for further migration (with probability,l-pysurvive)he dies). The migration thus implemented is a random walk with constant probability of stopping. We postulate our mites to search blindly and "bump into" orange locations. It is important to note that emigration is not continuous but occurs only at certain points in the local cycle. We have shown (Zeigler, 1977b) that continuous iration is unlikely to stabilize a locally unstable system such as Huffaker's.

15. Predator migration is carried out exactly as prey migration with respective parameters pdrem, pdsurvivc, p(ld:i f, pddi with the followinlg exceptions: 1') Predator migration is initiated when a local maximum in predator density is reached 5') After time Ts has elapsed, the individual is added to the popula.ti.o of his assigned cell unless the prey population is below eqprey, the equilibrium prey level computed from Eqs. lb) and lc). If t}he prey population is below eQpreyf a small-population stochastic iulocel takes effect. With probability drem', the invading predato.(-: remaincs otherwise it is migrated as in step 3). A predator that remains dies with probability d". If the predator lives, it creates another predator if there are at least c" prey and cx are used up as a result (where X is the current prey population size). The foregoing decision sequence is instantaneously computed and with TS = 1, the predator(s) are returned to step 5').

Note that although the small-population submcodel is a stochastic version of the deterministic Lotka —Volterra mnodlel. usecd for large numbers, it may be a crude summary of the local interaction in these circumstances, thus optimal settings of the primed parameters may bear little relation to their deterministic counterparts. Indeed, our methodology suggests that a second level spatial characte.:iz — ation of each orange could be built. Such a model would be tuned to mo.rcA fTn.e]y structured local observations and a sin plified version would replace (or perhaps turn out to be identical with) o current base submodel. In our general term — inology, no such experimental frame currently belongs to ~<., though one might: be forced to create such a frame, if the current models (in which the smallpopulation submodel participates) prove unable to match the data gathered within existing frames. (See Appendix.) 5.2. Discrete Event Lumped Model Our first lumped model (illustrated in Figure 3) is a discrete event ve-:sion of our base model. The migration is unchanged but the local interaction is escribed in summary transition function form obviating the necessity for step-by-step simulation of the differential equations. For example, consider the food-prey submodel, Eqs la) and lb) with y = 0. For positive r(O), it is possible to solve analytically to find r(t) = r(O) - u (x(t)-x(O))..2a) a x(t) = x(0) eat...2b) where a = b-d (nrt prey growth rate). The time for r to reach O is given from Eq. 2 by

!PREYI [REY PREY predator migration M~.AXPD / I X VAP PREDI ~ _/ / { \ A MAXPYT... O INTDECT ER GRT PREY3 PYDECT RE PREY food depleton RENEWT | food renewalI x[ prey population sequential y predator population tstate state r food amount e elapsed time Figure 3

17. ar(O) 2 ux(O) and the prey population at that time is: x(T) = ar(O) + x(O).. 3b) The discrete event model keeps track of the values r,x, and y for each cell. If at some time t, a prey individual migrates to an empty cell with food znounl Eq. 3a): with x(O) land r(O) = r is used to schedule the subseequent em:igr-ati occur at t ~ T). When the enigration event occurs, Eq. 3b is used to update the prey population, an6 of course, the food amount is set to z-cro. The consequent prey die-out is scheduled to occur inJ elapsed time (i/dl-iM)4[[..)p At every subsequent immigration the cell state is updated. Suppose that. a time e has elapsed since the last immigration. Then using Eq. 2 with t = e, and r(O), x(O) being the values at the last update, the correct state pertaining just before the immigration is computed. To the prey number so computed we add 1 to account for the immigrating individual and then use Eq. 3a) to reschedule the emigration event. It can be showm that this discrete event algorithm exactly reproduces the behavior of the base model prey-food interaction. The addition of the predator is handled in principlein the same way, except that the scheduling and update functions cannot be obtained analytically but can be approximated with a once-and-for-all sampling of the trajectories generated by Eq. 1 (Zeigler, 1977c). We shall provide a brief description of this process. Consider the case where prey have colonized an orange but have not yet exhi the food. Then Eqs lb and lc) reduce to the Lotka-Volterra dynamics and one easily obtains Jhe equilibrium prey and predator isoclines, namely

predator eqprey (eqprey, y3 ) /Th() iood >0 (x0gd,eqpredr) 5 eqpred (x~,e pred) T I / ~'........ 60 prey predator eqprey food = O (x,0) 60 prey Figure 4

18. x* = eqprey - d'/c' y* = eqpred - c/a These isoclines divide the plane into four regions as shown in Figure 4a).. A typical trajectory initiated by a predator immigration is segmented into an initial joint growth phase T(1) ' a predator growth phase T(2) and a joint decline phase T(3) Generalizing Huffaker's observation to predatorso ve postulate that predators emigrate at mnaximum predator population at the end of the T 2) phase on the eqprey isocline and at: the end of the -(3) pha- e whnl' the prey minimum is reached (ei`::r a cro:.-Ssh occurs:if pyredators. are nlumc-rouc e.no or the eqpred isocline is reached). Any prey left are assumed to remain on orange and subsequently take part in the standard food-prey interaction. Scheduling and updating for each of these boundary crossings was done by use of tables generated from a CSM 1 simulation of Eq. 1 and shown in FiguLre 5 These tables are interesting in themselves; to our knowkedge, they repriesent the first such global study of Lotka Volterra dynamics. The parameter T shown is the period of the cycle obtained by linearization around the equilibrium point, T = 2r//ad. The time to cross from one boundary to the next is approximately T/4 near the equilibrium but declines rapidly as inLiial populations increase. A disadvantage of the generating tables by simulation is that it must be done potentially anew for each set of parameters. This makes it important to be able identify the local interaction parameters before all others (Section 7.1). The theoretical basis for discrete event representation of systems is given by Zeig er (1977c).

(maximum prey) eqprey 300 yC(]) mnumber nf invading predators 60 0 60 300 xo (prey at time of predator invasion ) (time to reach eqprey maximum prey) ____ ~~~T/4 05 6030=5 0~~~~X

eqprey (maximum predators) I fooiti::> C t <-~~food>, ()X 10- 1/./ —food I0 0 100 500 x(D (moimrnum prey) YO B eqprey (time to reach maximum predators) __ T/. 5T food> 0 food=O O I00 500 xg (maximum prey) Figure 5

x' eqpred Ye (prey left after d (predator left predation) after prey anihilated) 60 eqprey food 0 food >0 food 0 I0 0 (mximm pretor) Figure y(5 (maxim'um predators )Figure

Y 2 eqpred ( time to reach of cycle) JOINTDECT 10d food =0 food > 0 0 -0 IF0 5O YF Figure 5 f

19. 5.3. Occupancy Models The overall occupancy model is; described in Pigure 6. As can be seen, the local description is reduced to a small number of states (empty, patch occupied by prey only, etc.) and the local dynamics are reduced to timed transitions from state to state. The migration component is also simplified by specifying neighborhoods for each species (for each cell) and migration is effected by means of independent Bernouli trials governed by specified probabilities and conditions at the cells in the neighborhood of a migration:active cell. (See Zeigler,1977b for a full explanation.) Cell spaces up to quite large sizes (we have commonly investigated 30 x 30 arrays [900 cells]) can readily be simulated in discrete event languages such as SIMSCRIPT. We are able to study by this means the spatial patterns which are associated with persistence and extinction as they are governed by the geometry of the space, the characteristics of the neighborhoods, and the settings of the other parameters. 5.4. Random Phase-Space (RPS) Models The RPS hypothesis assumes that the cells in a given state are uniformly distributed in both space and phase (elapsed time in the state) at all times. On the basis of this hypothesis we may derive the differential equation system shown in Table 6. These equations are simple enough to be solved for equilibrium isoclines and thus give qualitative information about how persistence is governed by the various parameters. Here persistence is judged relative to prespecified extinction levels such that if the prey and predator occupied cell fractions fall below these levels the system is

GRT2 '- DEC~2 0DECT1 E o t p Ih R cell E. Wait 0. Wait 1. Hold (GRT1) Y cell' in N1(cell), cell' in 0 -+ cell' in 1 Pl1 Hold (DECT1) go to E 2. Hold (GRT2) V cell' in N2(cell), cell' in 1 -+ cell' in 2 P2 Hold (DECT2) go to 0 R. Hold (RT) cell in E -> cell in 0 Figure 6 go to R

Table 6 Random Phase - Space Model Let x = fraction of cells in state 1 (somie prey) y = fraction of cells in state 2 (soeipe prey, some predator) z = fraction of cell.s in state 0 (empty, soome food) u = fraction of cells in state F(einpty, no food) Then dx x PiN1 P2N2 d-t T. + zx T yx T dy = -T+ yx P2N2 dtVT2T ~u -1-(x-l-y+Fz) dtT T T2 2 22 dz u P +1 y dt RT T1 T2 where T1 = GRT1 + DECT1, T2 = GRT2 + DECT2 Equilibrium Relations Equilibrium fraction nR predator (food,prey/pred) some predator (foo d prey, pred) x*(avg. prey cell) 1-(p1N1) /l+RT T (P y*(avg. pred cell) 0 l-(pN1) -(P2N2 (i4 RTT - ) l+p2N2T2 (PNl1) T1 z*(avg. food cell) ( 1) (P1N1 +y* P2N2T2 (pN) T1 u*(avg. utilization) l-(x*+z*) l-(x*+y*+z*)

20. assumed to go extinct. Moreover, the equations can be easily simulated to generate the associated dynamic behavior. The predictions thus made can be matched against the behavior of the simulated occupancy models, to the advantage of both model types (see Validation, Sec. 7). 6. Organization of Models We display the hierarchy of models in Figure 7. Models are organized in a manner parallel to that of experimental frames. Model M is on the same plane with model M' if M is a subcomponent of M '. Models of greater degrees of abstraction (lower level planes) are derived from more refined models by simplification procedures based on aggregation, coarsening and discreteeventization (as discussed in Section 3). The criterion which justifies the simplification is that the mappings involved be homomorphisms. We have given extensive expositions of such model relations (Zeigler, 1976a) and their use in organizing models (Zeigler, 1977b). For an exposition and example in the compartmental ecosystem context see (Zeigler 1976b).

Elocal,total \ <,\ I',. food(Ipreypred food prey predator local. - replenish- migration migration loca vlacWL-z ment. f ood 5,prc, /Jffc d Eglobal LUMPED global,tota od, preypred food,prey occupancy 1d E eS X \food OCCUPANCY Eoccupancy, total ~ ~ -__- \ RANDOM HASE-SPACE Figure 7

6.1 Relation of Parameters Basically a homomorphism betwiee-l models is; a corresponridencee between 1..-a'ilr state spaces in which corresponding states transit to corresponding states and yield corresponding outputs. Such a relation usually implies a correspondence between parameter values as well, so that parameter settings of a lumped model may be completely determined by those of a more refined morphic preimage. Indeed, we have given an ex-: mrnple of such a parameter coriespondencce jrl our derivation of the discrete event lumped model (Section 5.2). The schedu.ling and update tables of the discrete event model are parameters - one can treat then as entities to be arbitrarily adjusted until the desired behavior is achieved. C the other hand, the morphism, which underlies the construction of the lumpedlmode uniquely prescribes these tables for each setting of the local interaction parame (b,d,u,d'tcc') The same concept of parameter correspondence is illustrated in the parameter complexes appearing in the RPS equations, as determined by the parameters of the occupancy model. (Table 6). To complete the chain, we should indicate how the parameters of the occupancy model are related to those of the more refined discrete event lumped model. In this case, however, a morphism cannot be established to hold strictly over the complete state spaces of the two models and we must be content with estimating average parameter value settings, or at least, ranges to which they can be bounded. We now outline how this may be done. In Table 7 we provide ranges for the patch life cycle parameters GRTi,DCT i = 1,2. The derivation is straightforward form Section 5.2. Although the

Table / Estimation of Occupancy Model Parameters Range of Values in Terms of Parameter Lumped Discrete Event Model Parameter Values ar0 GRT1 (Growth Time of Prey Colony) 0[, a un 1a ar DCT1 (Decay Time of Prey Colony) [ 0 d u 2a GRT2 (Time to Maximum Predator [0 o,T/2 ] (T = - ) Population) Iad DCT2 (Time to extinction of Predator [ 0,T/4 ] Population measured from Maximum Population Point) From Fig. 8 with number of samples ar0 p1N1 (Effective Prey Colonization = (1 - pyrem) pysurvive and u -- Neighborhood) random walk parameters means pydif, ar From Fig. 8 with number of samples c ar0 p2N2 (Effective Predator.- (1 - pdrem) * pdsurvive Colonization Neighborhood) c u and random walk parameters meanpdsearch, pdd pddif2

22. migration mechanisms in the lumped and occupancy models are not directly coinpar they can be matched through the notion of effective neighborhood.. The effectiv neighborhood of a species is the expected number of cells colonized in a migrat episode, given that all cells in the space are colonizable. In the occupancy m this is just piNi for species i. In the more refined model, the effect:ive neighborhood is the expected numiber of distinct cells accessed in a migraticon episode. Figure 8 plots the number of distinct cells accessed versus the numlbt of samples taken from the random walk distribution (with parameters typical ill I case of extended coexistence). ";'he maximum number of distinct cells acce.sc;ed cf be estimated by noting that the random walk distribution with mean search time and diffusivity d in one dimension appears (from simulation) to be normally distributed with standard deviation a = d/f-. Thus, the number of distinct ce S rises at first in proportion to the number of samples; then it approximates the number contained within a radius 3a of the active cell as the number of sam-lple increases to moderate values. (Theoretically it continues to rise very slowly be: this point.) Now the number of samples in the migration episode is just the numbi migrants and can be bounded above as indicated in Table 7. Combining this numbei with Figure 8 yields the effective neighborhood bound. Finally, we note that yet more refined models can be postulated which would place constraints on the parameters of the base model. Thus in our base model, the prey death rate in a patch d, the probability of survival, pysurvive, and the search time parameter, meanpysearch, are independently adjustable. If we postulate that a migrant survives only if his lifetime

n X C50. 40 ' 30 r meonpy search I6days a) C20 p, t pydif ID.: E Z 1I O0 I _1 1....i. 1 ~~ L....L__ _ 0 100 200 300 400 500 600 1000 HpdH4 Number of saumples p.i-Figure 8y Figure 8

23. exceeds his search time, where these random varJ ables are independent and exponentially distributed we derive the relation: 1 pysurvive d _ _ a pyuri_ 1 + domean pysearch Parameter values obtained after adjustment can be checked against this relation. Large discrepancy might indicate dependence of the variables or give cause to reexamine the model structure and/or parameter settings.

24. 7. Applicability of Frames to Models A frame E is applicable to a model 1M zif the conmpare variables specified by E are 1) included among the descriptive variables of M or 2) may be obtained from the descriptive variables of M by aggregation or coarsening. Figure 7 depicts the "core" of the appticability relation (see also TabJI.co 1 and 5). "Core" is used here because one can infer applicability to higher plane models of frames applicable to lower plane models (Axiom 8; Zeigler, 1977b). Roughly, if E is applicable to M, then M can potentially answer the questions of interest in E. M is valid for 4t in E if M can reproduce Xt (E), the data collected up to time t in E. Viewed another way, E applicable to M means that the data D (E) may be employed to identify the parameters of M, i.e., the parameters may be adjusted until a best fit with the data jj (E) is obtained. 7.1. Parameter Identification As shown in Table 5 the experimental frame and model organizations made it possible to identify the parameters in a sequential manner, thus greatly reducing the search space at each stage. The parameters relating to: local foodprey interaction, local\food-prey-predator interaction, prey migration, predator migration and finally predator-prey (small-population) interactionwere adjusted in this order. The test of such a procedure is that reasonable fits to the data are obtainable at later stages by holding fixed the parameters identified at earlier stages. When acceptable agreement at later stages cannot be obtained, this may indicate that the prerequisite independence assumed for earlier adjusted parameters does not hold. In terms of experimental frames, the control conditions of a frame may not in fact hold. Indeed, it often

is implicitly assumed by modellers that certain global interactions can be ignored in certain circumstances, and this may turn out be unjustifiable. In our case, the E total frames assume that migration effects have been nullified, the justification for which lies in Huffaker's verbal account of migration episodes accompanying the time series data. If there is reason to doubt that the control conditions of an experimental frame are not satisfied, a readjustment of parameters may be attempted.. To the extent that such a readjustment is Sm.1l, the': derormposi.lion into cxperimentrTl frames will have been beneficial. In the Appendix, we report on the parameters identified in some of the key experiments and cross compare models in this regard.

26. 8. Summar The concepts discussed in the paper are summarized as follows: ' ',~ t' ~is the triple of conceivable experimental frames, rea:. system data and possible models, respectively, underlying the modelling study of the Huffaker universes. <&t =1961 is the subset of frames for which data had been collected until 1961. t = 1961(E) denotes the data collected within frame E until 1961. t = 1977 is the subset of models considered until 1977, the current time by the present modeller. An experimental frame E in C specifics a (compare,control) variable pair. A data element of a frame E is a time series of compare variable values obtained under conditions where control variables are kept at zero levels. Frames are partially ordered by the derivability relation; E < E' means that data elements of E are de-:civable from those of E' by employing selection, aggregation and coarsening operations. Each frame may be assigned a degree of abstraction equal to the minimum number of operations required to derive it from a fixed most inclusive frame. Models may similarly be partially ordered by use of morphism relations. A homomorphism is a mapping from a refined model to a coarse one which preserves the transition structures. A homomorphism induces a mapping from the parameter assignments of the finer model to those of the coarse one. If a frame E is applicable to a model M, this means that the behavior generated by M can be interpreted as data within frame E. One interpretation of this fact is that the real system data collected within E can be employed to identify the parameters of M.

27. 9. Discussion There are a number of levels at which the integrated approiA'h to modelling illustrated here may be discussed. We briefly consider some of them. 9.1 Large Scale Multi-Faceted System Modelling Our formalism has been constructed from a general starting point - the theory of systems and its specialization to modelling and simulation. Thus, it is aimed for applicatior to 'ltargqe scale" systems ins ge3ne ral. We have placed the large scale in quotation marks to signify our belief that "large scaleness" is a matter of approach rather than of fact. Indeed, a real system is called large scale precisely when one recognizes that to deal with it successfully requires the consideration of many factors and aspects. There are some systems which strikingly have this characteristic - environmental systems, urban systems, etc., that are indeed large scale. But "micro scale" systems such as the biological cell are equally comiplex, when examined in all their facets. Thus we propose the term "multi-faceted" to connote the systems (viewpoint) we are adressing. In this paper, we have illustrated our approach in a particular ecosystem context. But some general points clearly emerge. These are: Simpler models can give qualitative and sometimes quantitatively accurate predictions. The RPS model gives good estimates of average cell occupancy fractions when its underlying conditions hold. More generally it may give correct

28. qualitative relationships (effect of paraxileter s(fttinlg;s) even when its quantitative predictions are inaccurate. Simpler models can be employed to check more complex ones. If the correspondence between models is known, behaviors of the models may be compared. This can be employed at: a) the development stage; if the simpler model is known to be correctly implemented (or does not require simulation), then the logic of the more complex model can be verified by comparison of model behaviors (this is an important special case of redundancy use for program verification, Bosworthl-, 1976). b) the prediction stage; the more the predictions of various models agree, the greater may be the confidence in the predictions. Where serious disagreements occur, confidence considerations may determine the choice of which to believe, or lead to the conclusion that more development is necessary. Complex models can be employed to validate simpler ones. Conversely, if the correspondence between models is known, a more refined model whose details are tied to a particular real system can be used to gain confidence in a more abstract but general model. Thus by validating our lumped discrete event model against Huffaker's data, and finding that our corresponding occupancy and RPS models produce matching behavior, we gain confidence in the abstractions employed to derive the simpler models, i.e., that patches, rather than individuals, are sufficient entities for analysis of persistence. Holling et al (1974) has employed a simulation model at the level of detail of our lumped discrete event model to check out the wider consequences of optimal control policies derived from a simpler analytic model.

29. Models may be introduced independently or derived from existing ones. It may be sometimes advantageous to construct a model from "phenomenological considerations" rather than fromn "first principles". However, erhe. a homomorphism can be established between a more refined model and such an ad hoc model, additional advantages of the kind indicated above accrue. In addition, if a base model is available on which to base construction of a lumped model, constructs may be suggested which would not have come to mind in a phenomenological approach (Whitehead, 1977). Needed experiments may be imrlied by the experimental frame organizationl The logical structure of the experimental frame organization may suggest conceivable frames in & that have not yet been realized to date (are not in ) ), and might not be thoughlt of in an unstructured experilf->ntal. approach. For example, data on the orange spoilage and prey-food interaction suggested by frames of the form E and Efood,prey/pred is foodlpreypred foodprey/pred missing and would be helpful to model construction and validation. Model and experimental frame or anzationls mlay be extended at all levels,. The multifaceted system approach explicitly recognizes that model construction and validation is a never-ending process. For example, as accuracy demands in some frame increase, it may be found that the current stock of models is inadequate to meet these demands. This may spur the formulation of new experimental frames, data acquisition within them and construction and validation of models which would guide the refinement of the original models so as to meet the increased accuracy requirements. This paradigm is illustrated in our finding that small-population interaction on a patch may play a more imporitant role in determining average population levels than was suspected originally. Development of a credible small-population sub-model could be based on a spatial model of the predator-prey interaction

30. on a patch developed from data acquired in an appropriately defined experimental frame. Complexity conetraints would prohibit incorporating such a spatial model directly into our local interaction moldel anld thlus simplifications would be sought perhaps resulting in refinements of the classical Lotka-Volterra model along lines developed by Hassel et al. (1976). An example of refinement at the other extreme of abstraction is givel by the incorporation by Gurney and Nisbet (1977) of fluctuation terms in the RPS model which enables it to predict equilibrium fluctuation magnittucldes from steady state population levels. 9.2 Ecosystem Modelling In this paper we have illustrated our large scale multi-faceted approach in a highly restricted ecosystem context. Having dealt only with two species and 3 trophic levels, we have only scratched the surface of the possibilities and problems that would arise in dealing with a realistic ecosystem. Yet extension of the experimental frames on the same plane of abstraction to many species would simply involve the specification of frames by pairs (A,B) where A is the

31. subset of species to be observed (whose descriptive variables are the compare variables), and B is the subset of species where influence is to be min:mizc-C (whose descriptive variables form control variables). A sublattice of frames represents the trophic structure of the ecosystem such that (A, B) is in the sublattice if, and only if, the A species are found at lower trophic levels than any of the 13 species. Competitive and cooperative stuctures may be similarly represented. In the same vein, we have hinted at more than one level of patch decompos Indeed, the spatial structure II f have a natural hierarchy, where patches isolated at one level of analysis are subordinated in larger patches at a highe level. Isolation of patches signified by frames bearing the "local" descriptor would then be possible at many levels, and the frame characterization would reflect this hierarchical structure. Finally, ordering of frames according to plane of abstraction is clearly extendable to many degrees of abstraction, representing many possible sirmp].ifi — cation and aggregation procedures. Aggregations in time may result in successi' alternate planes of differential equation and discrete event models. Aggregatib over trophic levels and/or over patch hierarchy levels are possible. Compartme] alization of species according to spatial and functional criteria is yet another source of abstraction-frame construction. Parallel to such an experimental frame structure would be the organizatioi of models intended to answer questions within various applicable frames.

32. 9.3 Modelling of Patch Structured Predator —Prey Universes We have shown in particular, how the organizations of frames and models look in the highly restricted universes of Huffaker. Yet these universes are rich enough to enable non trivi.al multilevel, multiforma.isrm model construction and validation. Thereby we have illustrated that there are advantages to integrated modelling, even in restricted contexts more or less amenable to conventional treatment. Some tentative conclusions concerning predator-prey coexistence ar e: Predators and prey can coexist indefinitely in patchy environments where, in homogenous environments the relation would quickly go extinct. The necessary characteristics of patchy environments are: a) largely isolated patches - population exchange is small but non negligible, b) patch life cycle insensitive to emigration and to immigration except at certain key points (prey colonization, predator take-over), c) large number of patches. According to our base model, in the Huffaker Universes large numbers migrate,so small population exchange can be achieved only by making migration hazardous; the effect of emigration on the native life cycle is small since emigration occurs only at certain key points; and the effect of immigration is small because of exponential prey growth after colonization and rapid extinction after predator take-over. Gurney and Nisbet (1977), suggest a model in which a) and b) hold, but

33. emigration is possible throughout the cycle. Despite this discordance with Huffaker's observations, their model can be made to fit the persistence cases as viewed in the Eoccupancytotal frame (means and standard derivations of occupancy counts). Indeed, their model assumes the random phase - space mode of operation, and under these conditions, their model equations are isomorphic with our RPS equations (modulo fluctuation terms and inessential extra states). Our cross comparison of models.(Appendix ") suggests that the RPS mode can notbe maintained with the number of patches employed by Huffaker, so that even though a RPS model can fit the data, it may not represent any base model which also does so. Two kinds of tests of the Gurney-Nesbit model are suggested in the spirit of the multimodel approach: a) construction and test of a base model satisfying their assumptions, b) conducting critical experiments to distinguish the alternative mechanisms. Acknowledgements The assistance of Ramon Guarldons in many phases of the present project is gratefully acknowledged.

REFE RENCES Bosworth, J.L., (In Press), "Models and Simulation of Programs: Application to Self Checking Programs", In: Frontiers in Modelling, Ed. by B.P. Zeigler. Gurney,W.S.C. and R.M. Nisbet, (In Press), "Fluctuations in Predator-Prey Populations in Patchy Environments". Hassel, M.P., J.H. Lawton and J.R. Beddington (1976), "The Components of Arthropod Predation", J. Anim. Ecol. 45, 135-186. Hassel, M.P. and R.M. May (1974), "Aggregation in Predators and Insect Parasites and its Effect on Stability", J. Anim. Ecol. 41, 567-594. Holling, C.S. et al (1974), "Project Status Report: Ecology and Environmental Project". IIASA SR-74-2-EC, Laxenburg, Austria. Huffaker, C.B. (1958), Experimental Studies in Predation: Dispersion Factors and Predator Prey Oscillations", Hilgardia, 27(14), 243-383. Huffaker, C.B., E.P. Shea and S.G. Herman (1963), "Experimental Studies on Predation: Complex Dispersion and Level of Food in an Acarine PredatorPrey Interaction", Hilgardia, 54(9), 305-329. Levin, S.A. (1976), "Spatial Patterning and the Structure of Ecological Communities", Some Mathematical Questions in Biology, VII, A.M.S. Maynard Smith, J. (1974), Models in Ecology, Cambridge University Press. Whitehead, B., (In Press), "Base-Lumped Neural Models For Human Pattern Recognition", In: Frontiers in Modelling, Ed. by B.P. Zeigler. Zeigler, B.P. (1976a), Theory of Nodelling and Simulation, Wiley, N.Y. Zeigler, B.P. (1966), "The Aggregation Problem", in Systems Analysis and Simulation in Ecolqgy, Vol. 4, ed. by B.C. Patten, Academic Press, 1976.

Zeigler, B.P. (1977a), "Structuring t:he Organizat:ion of Partial Models", Int. J. Gen. ys. (In press). Zeigler, B.P. (1977b), Persistence and Patchiness of P'rcedator-Prey Sy';tems Induced by Discrete Event Population Exchange Mechanisms", J. Thoer. Biol. Zeigler, B.P. (1977c), "Systems Simulateable By the Digital Computer: Di:scr.ete Event Representable Models", Logic of Computers Technical Report, University Of Michigan. Zeigler, B.P. and A.G. Barto (1977), "Alternative Formalisms for Bio --- and Eco-System Modelling", Simulation.

APPEN DIX Some Results on Estimated Parameter.; and Model Cross-Comprsoin It is beyond the scope of this paper to describe the complete results of simulation studies of the constructed models. We shall describe, however, the results relevant to the model reproduction of predator-prey coexistence as observed in the (one) 1958 universe which manifested this phenomenon. In general, all of our models are capable of qualitatively explaining the observed persistence, and our more refined models are able to give quite close quantitative agreement as well. We proceed to describe the results of the identification procedure outlined in Section 7.1. Refer also to Table 5. A.1 Local Food-Prey Interaction -1 — 5 The parameter settings b = 0.55 day, d - 0.33 day and u = 1.3 x 105 (fraction of orange surface per day) wereestimated by employing the initial part of data element 58, IB, Fig. 7 (in frame Ef docra/tro d). The quantities food,prey/pred initial prey size, initial orange surface exposed, maximum prey size, time elapsed to maximum prey were employed in Eq. 2 to estimate a(=b-d) and u; the slope of decay from maximum prey size was used to estimate d. At the above settings, the maximum prey size on an orange is u/a~ 2 x 104 mites per orange equivalent (so for example a 1/10 exposed orange area can support 2000 preys at the perigee of growth). A.2 Local Food-Predator-Prey Interaction The parameter settings d' = 0.30 day4, c = 0.05 day and c' = 0.006 day1 were estimated employing data element 58,IIA,Fig. 9 (in frame Ef~ d'd ad) Employing a CSMP simulation of Eq. la,b,c) we adjusted the parameters d', c and

c' so as to fit as closely as possible the prey and predator curves in Fig. 9 (Huffaker, 1958). With the estimated parameters, we have the equilibrium prey and predator levels as 50 and 5 mites on a patch, respectively. A.3 Prey Migration The data necessary for identifyring prey ming)ation parameters in the absence of predators is available only for the hazard free universes in the Huffaker 1958 study and the complex universe of the 1963 study, but not for the 1958 universe i n which p /- ypredator coexistence was achieved.?,hnployJng hazard free universe data sets 58, IB, Fig. 7 and 58, IC, Fig. 8 (in frarmev Efood,reyped) we adjusted the parameters pyrem, pysurvive, meanpysearch,4 and pydifl of our discrete event lumped model so as to have the SIMSCRIPT generated curves match the data curves as closely as possible In mnaximu.m prey produced and number of prey maxima produced in the experimental interval Estimates obtained were pyrem = 0.9, pysurvive = 0.9, meanpysearch = 0.1 days pydif1 = 20. Employing the data set 63 E-2 Fig. 2 for the prey-food inter — action in the complex 1963 universe, we estimated in the same manner that pyrem = 0.9, pysurvive = 0.5, meanpysearch = 13 days, pydif1 = 0.3 and ydif22 = 0.15. Thus as expected prey mites in the complex universe take rntch longer on the average (13 versus 0.1 days) to cover much less distance (/13 x 0.3/-: 1.0 versus /671 x 20~s 6.0, see Section 7) than they do in the hazard free cases. A.4 Predator Migration Predator migration parameters pdrem, pdsurvive, meanpdsearch, and pddif were adjusted in the discrete event lumped model so as to fit as closely as possible th data element 58, 11 I, Fig. 18 representing the 1958 universe in which coexistence was established. The settings of the prey migration

parameters were those determined from the complex 1963 universe just described, (Subsequent trials with deviations from these settings did not significantly improve the results.) The predator migration parameters were initially set equal to those of the prey and a fairly broad neighborhood of parameter assignments centered on the initial settings was investigated. It was found that coexistence is robust in this neighborhoud in that most simulation runs ended with both predators and preys still around. However, it did not seem possible to achieve very close quantitative agreement. We noticed that the predator occupied cell fraction was too small aud this seenmeTd to be due to the fact that in our original model, predators invrading patches of low prey density (less than eqprey) were always returned immediately for continued migration. It thus appeared that predator invasion of low density patches was a significant process and we accordingly modified our small.. population submodel to its current form. With this modification we were able to bring the statistics shown in Table A.1 generated by the simulation quite close to those of the data. Although the averages agree quite well, the model. overestimates the prey maximum considerably, which may point to a further needed modification. (In analogy with the predictions of a Lotka —Volterra model, the overshot could be the sensitive result of too low an initial predator population, and thus not an intrinsic model shortcoming.) The best fit parameter settings are indicated in Table A.1. In Table A.2, the same data is analyzed from the cell occupancy point of view (frame EoccuPancytotal) [Note the model in question is the lumped food,prey,pred discrete event model not the occupancy model; the occupancy states can be computed from the finer population count information. ]

Comparison of Data and Lumped Discrete Event Behavior in Frame global' total Frame Eypred [the Case of Predator —Prey foodprey,,pred Coexistence, 58 It I, Fig. 8] Density Data Modeli* maximum prey mi ss;ing 4600f (predators absent) average prey missing 2400 (predators absent) maximum prey 2000 3500 average prey 900 730 maximum predators 50 46 average predators 12 13 +all densities quoted in mites per orange equivalent (Huffaker, 1958). *Parameter assignments are: b = 0.55 day pyrem = 0.9 pdrem' = 0.3 d = 0.30 day' pysurvive = 0.35 d" = 0.0 u 1.3x 10 meanpysearch -- 13 days c" = 10.0 d'= 0.30 day pydif1 = 0.3 c = 0.0 -1 c = 0.05 day c'=.005 day pdrem 0.6 pdsurvive = 0.5 meanpdsearch = 14 days pddif1 = 0.2

Table A. 2 Comparison of Data and Lumped Discrete Event Model Behavior in Frame EoccuPaucy total food,prey,pred Cell occupancy Dat a Model average prey cell (state 1)+ 17 28 average pred. cell (state 2) 131 15 standard deviation/prey cell 16 19 standard deviation/pred. cell 11.1 a prey cell is a cell occupied by at least some prey but no predat:or a predator cell is a cell occupied by at least some predator measures the amplitude of oscillation considered as a fluctuation about the average (Gurney and Nesbit, 1977)

As can be seen, the statistics front mode]. and data are remarkably cl.ose, save for considerable overestimat:-ionl ill the average prey cell count. Thi s is understandable in view of the maximum prey population overestimation. It should be noted that the average occupancy counts are not necessarily correlated with the average population counts. As we have noted, the somewha independent occupancy perspective was useful. in diagnosing s shortcoming of the model. A.5 Occupancy and RPS Models Employing the parameter values of Table AdI, we can determine correspond parameter values for the Occupincy model, us-ing the relations of Table 7. In order to explore the behavior of the occupancy model in this space, we fixed all but the migration parameters at the extremes of their ranges and sampled the model behavior for allowable assignments of the latter parameters Employing the equilibrium relations in Table 6, we can uniquely determ:ilne the effective neighborhoods plN1 and P2N2 of the RPS model required to reproduce the occupancy averages of the data (58k II I, Fig. 8) shown in Table A.2. As shown in Table A.3, these are within but at the lower end of the ranges computed from Table 7. However, simulation of the occupancy model with these parameter settings resulted in quick elimination of the prey. Only when the effective prey neighborhood was considerably increased and the effective predator neighborhood considerably decreased was coexistence obtained in 10 x 10 cell array. (Halving the predator neigborhood was sufficient for coexistence in a 30 x 30 cell array. The 100 cell array is more representative of the 120 cell 1958 universe.) The effective nejghborhc obtained in the way are still within the ranges computed from Table 7. However, the occupancy averages obtained from the occupancy model for both predator ani prey in these cases tend to exceed those of the lumped discrete

Tabl]e A. 3 occup ancy total Cross-Comparison of Model Behavior in Frame Efocpcy, food,prey,pred Effective Prey Effective Predato-:i Occupancy R'S+ Neighborhood Neighborhood Model Modei avg. avg. avg. avg. prey pred. prey pred. P1N1 P2N2 cell ell cell cell. c[0,25] ~[0,12] 4 8.0 x 10 Cxt inct 17 17 30 x 30 extinct 4 4 10 x 10 extinct: 30 x 30 44 9 30 18 24 24 1.0 x 10 60 28 48 50 Data 17 11 17 11 Lumped Discrete Event 28 15 28 15 +Other parameter values: GRT1= 20 days, DCT1 = 20 days, GRT2 5 days, DCT2 = 2 days, RT = 44 days From Table A.2

event model and the real system data. In sum, this between model comparisort seefns to indicate that the random-phase condition is only approximlalelt.e'ly bI:i.ng satisfied in the ltwiped discrete event model and the real syst:em. While the occupancy and RP'S models predict that coexistence is possible within the allowed parameter space, they do not do very well in predl.icting the observed occupancy cell. averages unless the number of cells is considerably increased.

CAPTIONS Figure 1. A data element of a frame Elobal to is sown ill Fig. ib) food,preyfpred (redrawn from Fig. 8, Huffaker, 1958). The universe consi.sts of four oranges embedded in an array of rubber balls (oranges are the darkened circles in Fig. la). The hatched initial portion of Fig. lb) is the data element belonging to the frame Elocal where the subset food,prey/pred w referred to by the "local" designation is the set of oranges indicated in Fig. la). Figure 2. Experimental frames organized according to planes of abstraction. Degree of abstraction increases from top to bottom. Nodes represent frames and lines (implicitly directed from top to bottom) represent: the derivability relation. Figure 3. Discrete-event representation of the base model. The discrete states are: E(empty,no food), ER(empty,food replenished), PREY(prey colony established), PREY' (prey colony at maximum size), PRED(predators invaded), iA)XPD(predator colony at maximum size). Scheduling times are: GRT(growth time of prey colony), PYPECT(decay time of prey colony), MAXPYT(time to reach maximum prey size after predator invasion), PAXPDT(time to reach maximum predator population from maximum prey population) and JOINTDECT(time from maximum predator to end of cycle). Figure 4. Typical trajectories in the prey-predator (x,y) plane with food r > '0 (Fig. 4a) and r = 0 (Fig. 4b). Figure 5. State update and scheduling curves obtained by simulation of Eq. 1. Symbols shown are keyed to Fig. 4.

Figure 6. The occupancy model. Discrete states are: O(empty,food replenished), l(some prey), 2(some predator), E(empty,no food). Figure 7. The organization of models.t= 1977. Also shown are the experimental frames applicable to the various models. Figure 8. The cumulative number of cells hit versus the number of samples from the random walk distribution with parameters typical in the case of extended predator-prey persistence. The numbers pd and XJ indicate upper bounds on the numnbers of 1)redators and preys emigrat.:i1g in a migration episode as estimate t-,L in Table 7.