COEXISTENCE OF SPECIES IN A FLUCTUATING ENVIRONMENT by Gordon James F i t z p a t r i c k B.A., University of Alberta, 1965 B.Sc., University of Calgary, 1970 H. i . S c , University of B r i t i s h Columbia, 1972 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS OF THE DEGREE OF DOCTOR OF PHILOSOPHY i n the Department of E l e c t r i c a l Engineering He accept t h i s thesis as conforming to the reguired standard THE UNIVERSITY OF BRITISH COLUMBIA June, 1977 © Gordon James F i t z p a t r i c k In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t o f the r e q u i r e m e n t s f o r an advanced degree at the U n i v e r s i t y o f B r i t i s h Co lumb ia , I a g ree that the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s tudy . I f u r t h e r agree t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y purposes may be g r a n t e d by the Head o f my Department o r by h i s r e p r e s e n t a t i v e s . It i s u n d e r s t o o d that c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l not be a l l o w e d w i thout my w r i t t e n p e r m i s s i o n . Depa rtment The U n i v e r s i t y o f B r i t i s h Co lumbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 ABSTRACT A dynamic model i n which multiple consumers of a s i n g l e nutrient may coexist i n a f l u c t u a t i n g environment i s given. Only one consumer can p e r s i s t i n a fixed environment, but coexistence may be produced by e f f e c t s of a fl u c t u a t i n g environmental variable on nutrient u t i l i z a t i o n d i f f e r i n g between consumers. An approximate s o l u t i o n i s given for the non-autonomous Lotka-Volterra-Verhulst ordinary d i f f e r e n t i a l eguations of the model together with h e u r i s t i c s u f f i c i e n t conditions for construction of a persistent multispecies consumer community. Computational examples demonstrate persistence of an i d e a l i z e d example community f o r perio d i c and random environmental f l u c t u a t i o n . Two further examples demonstrate that environmental f l u c t u a t i o n can produce coexistence when environmental variables, standing crops, assi m i l a t i o n e f f i c i e n c i e s , primary productivity, u t i l i z a t i o n rates, and re s p i r a t i o n rates are comparable to a t r o p i c a l grassland, and an olig o t r o p h i c temperate lake. The s e n s i t i v i t y of model solut i o n s to fu n c t i o n a l v a r i a t i o n s of the component species may be ra p i d l y and accurately calculated. This allows the i d e n t i f i c a t i o n and estimation of unknown species f u n c t i o n a l responses from time s e r i e s data of biomasses and a measured environmental var i a b l e . Unknown functions of an environmental v a r i a b l e are approximated by a Tchebycheff polynomial expansion i n that i i i v a riable. Unknown c o e f f i c i e n t s of these expansions are the parameters of the model. These parameters are determined by the unconstrained minimization of the squared deviations of the logarithm of biomass observations and model d i f f e r e n t i a l eguation solu t i o n using a Quasi-Newton algorithm. This l e a s t sguares estimator was applied to a one year biomass time s e r i e s of four zooplankton grazers, phytoplankton, and average lake temperature of a small o l i g o t r o p h i c lake. Application of the model to t h i s grazer zooplankton community gives evidence of p a r t i a l s t a b i l i z a t i o n due to environmental f l u c t u a t i o n i n a natural community. I t i s concluded that environmental v a r i a t i o n , which i s often assumed on t h e o r e t i c a l grounds to be d e s t a b i l i z i n g , should rather be considered as one of the bases of community persistence. iv ACKNOWLEDGEMENT The author wishes to p a r t i c u l a r l y thank C.J. Walters for hi s continued encouragement, suggestions, and c r i t i c i s m , and Caroline F i t z p a t r i c k who has followed the ups and downs with me. Thanks are also due other members of the candidates committee, P.A. Larkin, and G.F. Schrack, f o r t h e i r advice and encouragement, and the research supervisor, A.C. Soudack, for his support and aid i n putting the thesis i n t o f i n a l form. This research was supported by an N. B.C. Scholarship to the author and by N.B.C. Grant number 67-3138 to A.C. Soudack. Unpublished data used i n t h i s t h e s i s were obtained i n a study supported by N.H.C. Grant number 67-3151 to W.E. N e i l l , T.G. Horthcote, and C.J. Walters of the I n s t i t u t e of Animal Resource Ecology. V TABLE OF CONTENTS Introduction .............................................. 1 Motivation ....«,...................<>.......... 1 Previous Work .......................................... 2 A New Concept .......................................... 6 Scope of Investigation ................................. 7 Chapter I A Model of Coexistence Maintained by Environmental Fluctuation .............................. 9 1.1 Single Nutrient Multiconsumer Varying Environment Model 12 1.2 General Bounds on Model Biomasses .................. 14 1.3 Deterministic I n s t a b i l i t y of Multiconsumer System 15 1.4 Approximations to a Persistent Multispecies Solution i n a Periodic Environment ........................... 15 1.4.1 Memoryless Nutrient Approximation ............. 15 1.4.2 Memoryless Nutrient Consumer System ........... 20 1.4.3 Bounds on Consumer Biomass Variation for Periodic Solutions with Periodic Environment ..... 21 1.4.4 Bounds on Biomass Means f o r a Periodic Memoryless Nutrient Consumer System .............. 24 1.4.5 Constructive Conditions f o r a Model with a Persistent Solution .............................. 28 1.4.6 S t a b i l i t y of a Related Time Invariant System ..29 1.5 Nutrient Growth, Environmental Fluctuation and Consumer Species D i v e r s i t y .......................... 30 1.6 Conclusions ........................................ 33 vi Chapter II Examples of Coexistence: Numerical Solutions ... 34 2.1 Example 1: Ideal 36 2.2 Example 2: Grassland 55 2.3 Example 3: Lake 61 2.4 Conclusions ........................................ 68 Chapter III Model F i t t i n g Procedure ....................... 71 3.1 I d e n t i f i c a t i o n Method 74 3.1.1 Environmental Function Modelling .............. 74 3.1.2 Continuous Environment Modelling .............. 76 3.1.3 Minimization Problem ................. 76 3.1.4 Solution of the Minimization Problem 77 3.1.5 I n i t i a l Parameter Guess and Obtaining a Good Starting Point ................................... 79 3.1.6 Termination of Minimization ................... 80 3.1.7 Evaluation of Parameter Estimates: Graphical Display .......................................... 81 3.2 I d e n t i f i a b i l i t y of Model: Noiseless Observations ... 82 3.2.1 Test Data 82 3.2.2 I n i t i a l Guess and Assumed Structure ........... 82 3.2.3 Minimization Performance ...................... 85 3.3 Estimation of Model from Noisy Biomass Observations 85 3.4 Estimation of Model from Slowly Sampled Observations: Indeterminate Parameters .............. 94 3.5 Estimation of Model from Slow and Noisy Observations: Given r ( v ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 3.6 Results ..102 v i i 3.7 Discussion ............... .......................... 103 3.8 Conclusions ..104 Chapter IV F i t t i n g Model to F i e l d Data 105 4.1 Data ...... 108 4.2 F i t t e d Sub-structures ........112 4.3 Comparison of Functions to Independent Measurements 125 4.4 Interpretation of F i t s ...125 4.5 Predictions of the F i t t e d Model 130 4.6 Conclusions ........................................ 135 Chapter V Conclusion ...................................... 139 Results of Investigation ...............................139 Suggestions for Further Work ........................... 140 Conclusions ............................................143 References .................................................145 Appendix I General Bounds on States ....................... 149 i . 1 Non-negative Biomasses ............................ 149 1.2 Upper Bounded Nutrient 149 1.3 Upper Bounded Total Biomass ....................... 150 1.4 Positive Biomasses f o r F i n i t e t ................... 153 Appendix II Deterministic I n s t a b i l i t y of Multiconsumer v i i i LIST OF FIGDRES-Figure 2.1 38 Figure 2.2 40 Figure 2.3 - .......44 Figure 2.4 45 Figure 2.5 .... 46 Figure 2.6 ... ....48 Figure 2.7 50 Figure 2.8 ..........51 Figure 2.9 53 Figure 2.10 ...54 Figure 2.11 ......57 Figure 2.12 . ..............58 Figure 2.13 . ... 58 Figure 2.14 60 Figure 2.15 63 Figure 2.16 65 Figure 2.17 .......66 Figure 2.18 67 Figure 3.1 A . ........83 Figure 3.1B 84 Figure 3.2 - ..........87 Figure 3.3A ..89 Figure 3.3B .-.-90 Figure 3.4 -.93 Figure 3.5A 96 iz Figure 3.5B 97 Figure 3. 6 A ..98 Figure 3.6B 99 Figure 3.7A 100 Figure 3.7B -.101 Figure 4.1 .109 Figure 4.2 -- .-111 Figure 4.3A ......................114 Figure 4.3B 115 Figure 4.4A 116 Figure 4.4B ... ....117 Figure 4.5A ....................... - ...118 Figure 4.5B - 119 Figure 4.6A .....................................120 Figure 4.6B .. ....122 Figure 4.7A ....... . - 123 Figure 4.7B ..........124 Figure 4.8 ---131 Figure 4.9 -.-134 X LIST OF TABLES Table I : Fixed Environments And Resulting E q u i l i b r i a With Time Constants 42 Table I I : Optimizer S t a t i s t i c s : Exact Observations ........86 Table I I I : Optimizer S t a t i s t i c s : Noisy Observations .......88 Table IV: Optimizer S t a t i s t i c s : Slow Observations .........94 Table V: Optimizer S t a t i s t i c s : Slow Observations, Fixed r(v) 95 Table VI: Optimizer S t a t i s t i c s : Slow Noisy Observations, Fixed r(v) 102 Table VII: Summary Of Substructure F i t s ...................112 Table VIII: Comparison Of Estimated Functions To Measured Rates ..................................................126 Table IX: Solutions Of F i t 9 with Constant v .............. 133 1 INTRODUCTION Motivation Ecology has been defined as the study of the d i s t r i b u t i o n and abundance of organisms. I t may be distinguished from natural history by the formulation and t e s t i n g of hypotheses about the functional determinants of natural patterns. Many techniques and models from other d i s c i p l i n e s have been applied to e c o l o g i c a l s i t u a t i o n s i n the hope that the p o t e n t i a l power of a mathematical deductive theory to unify perspective and reveal unsuspected r e l a t i o n s w i l l be as successful i n ecology as i n the physical sciences. I t appears that i f e c o l o g i c a l systems have d i s t i n c t i v e and important r e g u l a r i t i e s , the forms of these have not been described or modelled i n any other science. These d i s t i n c t i v e elements must therefore be recognized and a mathematical theory b u i l t which encompasses them at i t s base f o r a usable e c o l o g i c a l theory to r e s u l t . Unfortunately the existence and detection of such general laws remains an a r t i c l e of f a i t h since there i s , at present, not one e c o l o g i c a l p r i n c i p l e without many serious exceptions i f the p r i n c i p l e i s formulated i n a testable manner. Nevertheless the application of what ecology i s known, guessed at, and imagined has become necessary i n order that the e f f e c t of d i f f e r e n t p o l i c i e s f o r management of renewable b i o l o g i c a l resources may be predicted, without a knowledge of fundamental e c o l o g i c a l processes, e c o l o g i c a l manipulations run the r i s k of attempting the impossible i n the same sense that 2 attempted construction of a perpetual motion machine i s poi n t l e s s i f the laws of thermodynamics are true. The lavs of thermodynamics are not true because they are reasonable and consistent, however. They are true because they make predictions which have been often confirmed by measurement and no exceptions have been found. Therefore a mathematical e c o l o g i c a l formulation must not only be reasonable and consistent, i t must also be testable. This i n v e s t i g a t i o n sought t o devise a dynamic multispecies mathematical model of simple but novel form which accounted for v a r i a b i l i t y and coexistence of species and could be tested against and estimated from biomass time s e r i e s data. Previous Work One may d i s t i n g u i s h between three broad categories of dynamic models: multispecies t h e o r e t i c a l models of which the o r i g i n a l Lotka-Volterra[ 1 ][ 2][3 J[ 4 3 equations are the h i s t o r i c a l and conceptual base, s i n g l e species estimable management models such as the dynamic pool models of Beverton and Holt[5] and surplus production models[6] which may be considered to be developments of the Verhulst l o g i s t i c growth law[7], and f i n a l l y the large and diverse set of computer simulation models used variously as hypothesis formulations, p o l i c y gaming t o o l s , concept f o c i and predictors. Computer simulation models of e c o l o g i c a l systems are usually constructed from a combination of experimental and observational knowledge of the components and t h e o r e t i c a l l y 3 guided i n t u i t i v e estimates of unmeasured i n t e r a c t i o n s . The problems of v a l i d a t i n g such a model a r i s e p a r t i a l l y from the i a r g e number of parameters i n comparison to the data a v a i l a b l e f o r v a l i d a t i o n , but more e s s e n t i a l l y from the f a c t that parameters may not be definable from the v a l i d a t i o n data. That i s , assuming a deterministic simulation model whose inputs and outputs may be observed without e r r o r , i t i s not usually possible to reconstruct the model from i t s inputs and outputs. Therefore the sense i n which a model can be considered validated by comparison to observations or even by acceptable predictions i s questionable. There may be an i n f i n i t e number of choices of parameters of the model which would produce the same model outputs. Consequently detailed simulation models such as r 8 9 ] were not considered since they are not constructed i n a manner which allows estimation from system behavior. Single-species models form the bulk of management models i n f i s h e r y and w i l d l i f e management. Neglecting i n t e r - s p e c i e s i n t e r a c t i o n s allows the construction of models with d e t a i l e d age-structure d i s t i n c t i o n s which can be i d e n t i f i e d and to which estimation procedures may be applied. V a r i a b i l i t y of populations may be accounted f o r by s t o c h a s t i c c o e f f i c i e n t s and estimation i s possible. They are however si n g l e species models and were not considered as an appropriate basis f o r multispecies modelling since the number of i n t e r a c t i o n s possible between age structured populations was f e l t t o be too vast to be estimated from observations. The development of mathematical multispecies models and theory from the o r i g i n a l work of Lotka and Volterra has been extensive and has led to many reexaminations of e c o l o g i c a l theory, p a r t i c u l a r l y i n r e l a t i n g productivity, d i v e r s i t y and s t a b i l i t y . Rosensveigf10] related productivity changes to s t a b i l i t y of predator-prey systems. The role of environmental disturbance i n l i m i t i n g niche overlap on a resource gradient has been illuminated by the work of Hay and HacarthurT11 ]# and many ideas have been c l a r i f i e d by the study of d i f f e r e n t i a l equations of the Lotka-Volterra-Verhulst(L-V-V) form. These models have been reviewed by May[1U]. These mathematical models have not been conspicuously successful as a d i r e c t source of testable general hypotheses and predictions however. The two main c r i t i c i s m s of these models for the stated goal are the treatment of population v a r i a t i o n s and the d i f f i c u l t y of estimating such models from f i e l d observations to give sensible f i t s and predictions. Regular v a r i a t i o n s i n populations have been treated as the r e s u l t of a self-sustained o s c i l l a t i o n or l i m i t cycle produced by a time i n v a r i a n t set of d i f f e r e n t i a l or difference eguations. These models have not been successful i n estimation from f i e l d observations[12 ], but have had some success i n explaining laboratory competition experiments^13 J and as guides to construction of simulation models./ Random va r i a t i o n s i n multispecies models have been recently considered[15][16] as the consequence of random perturbations to 5 the c o e f f i c i e n t s of a d e t e r m i n i s t i c a l l y stable set of L-V-V eguations. In these treatments the coexistence of species i s considered to be l i m i t e d by the variance of the stochastic s o l u t i o n about the deterministic mean s o l u t i o n , species becoming extinct due to an expanding variance of s o l u t i o n . The estimation from observations has no developed methods and the properties of such equations are s t i l l under inves t i g a t i o n . The assumption of uncorrelated random disturbance as a model of the e f f e c t s of environmental f l u c t u a t i o n may be questioned. Species assumed to be i n competition are i n part i n the same environment and competinq f o r the same resources so vould be expected to be affected i n a correlated manner. The approach of Long[17] and Long et al.[18] of replacing parameters i n the l o g i s t i c equation by deterministic functions of environmental variables and including some in t e r a c t i o n between species meets some of the c r i t e r i a of the stated goal. In [17] a herbivore-plant-nutrient model was given where carrying c a p a c i t i e s , K, for each species were functions of food density. I n t r i n s i c rate of growth, r, was a function of temperature and K. This model was f i t t e d with some free parameters to part of a time s e r i e s of Daphnia qalatea numbers, chlorophyll "a1* concentration, n i t r a t e concentration, and water temperatures, and the model used to predict the remainder of the time s e r i e s . The modelling of species at the same troph i c l e v e l [18] did not, however, include any provision for competition and the s p e c i f i c assumptions of functional forms i n t h i s approach 6 r e s t r i c t s the generality of the model. The response stucture -simulation approach of Haguire[19] permits modelling of the form of environmental e f f e c t s with great generality, but does not lend i t s e l f , by the graphical nature of synthesis of environmental response surfaces, to estimation from limited observations by a computational process. The format of model f i n a l l y envisaqed was to be a set of L-¥-V equations with c o e f f i c i e n t s which are environmentally dependent functions of measureable environmental variables. Yet t h i s modelling approach leads i n general to a parameter estimation problem which i s p r a c t i c a l l y i f not t h e o r e t i c a l l y i n t r a c t a b l e . A Hew Concept The requirement that a set of L-V-V eguations be stable when c o e f f i c i e n t s are not time-varying reguires that the i n t e r a c t i o n matrix be of f u l l rank. A l l previous modelling approaches to a system with time-varying c o e f f i c i e n t s have assumed that to obtain a s t a b l e time-varying system, the corresponding fi x e d or deterministic system must also be stable. From a phi l o s o p h i c a l point of view May fs approachf151 considers environmental f l u c t u a t i o n s as an a l i e n , that i s , independent of the system, disturbance to a stable deterministic system. Yet i t seems just as reasonable that a community p e r s i s t i n g i n a f l u c t u a t i n g environment w i l l have adapted to that environment and may actually require the environmental disturbance to p e r s i s t . To t h i s end the construction of a model 7 which was d e t e r m i n i s t i c a l l y unstable, but stable i n a f l u c t u a t i n g environment was attempted. Removing the condition of deterministic s t a b i l i t y allows the instantaneous i n t e r a c t i o n matrix to be singular, thus reducing the number of free parameters necessary to characterize the i n t e r a c t i o n matrix, reducing and s i m p l i f y i n g the estimation problem to a p o t e n t i a l l y p r a c t i c a l s c a l e . Preliminary numerical experimentation established that a renewable resource p a r t i t i o n i n g process could be devised which would allow persistence of multiple f r a c t i o n s i n a f l u c t u a t i n g environment and only one i n a f i x e d environment. The i n v e s t i g a t i o n then concentrated on putting t h i s process into an e c o l o g i c a l l y interpretable model. Scope Of Investigation Chapter one describes the structure of an elemental model which i s d e t e r m i n i s t i c a l l y unstable, but can p e r s i s t i n a f l u c t u a t i n g environment. Analysis of t h i s model gives an approximate s o l u t i o n with h e u r i s t i c s u f f i c i e n t conditions for the s o l u t i o n to e x i s t , and constructive conditions f o r an environmental f l u c t u a t i o n s t a b i l i z e d community. In Chapter two numerical solutions were used to investigate some of the persistence properties of the model f o r an i d e a l i z e d system with both per i o d i c and random environments. Further computational examples e s t a b l i s h that persistence can be obtained at growth rates and standing crops comparable to two natural systems under a periodic environmental v a r i a t i o n 8 comparable to the respective natural environmental v a r i a t i o n . Chapter three describes an i d e n t i f i a b l e parametric form f o r unknown general functions of an environmental variable which comprise the parameters of the model. The i d e n t i f i a b i l i t y of the model from noise free data i s v e r i f i e d using a Quasi-Newton minimization procedure to obtain a least squares f i t of model d i f f e r e n t i a l equation s o l u t i o n to a r t i f i c i a l sampled observations. The performance of the minimization procedure i s described i n some d e t a i l to give guidance to other a p p l i c a t i o n s . The effectiveness of the technique as an estimator of parameters when observation errors are present i s tested on a r t i f i c i a l sampled observations. Chapter four describes the application of the estimation technique and model to f i e l d data from the grazer zooplankton community of a small lake., The model and estimation technique allow conclusions to be drawn about that system. Suggestions f o r extensions of the model, use of the concept of a fl u c t u a t i o n s t a b i l i z e d system and uses of the model and extensions are given. CHAPTER I A BODEL OF COEXISTENCE MAINTAINED BY ENVIRONMENTAL FLOCTDATION Temporal v a r i a t i o n s i n environmental variables such as r a i n f a l l , l i g h t i n t e n s i t y , and temperature are ubiquitous i n natural ecosystems. Events i n the l i f e of many organisms are correlated with the more predictable changes such as da i l y and seasonal v a r i a t i o n , and l e s s predictable vari a t i o n s from these average patterns, such as an unusually wet summer, favor some species and suppress others. A community p e r s i s t i n g i n a time varying environment i s affected by that v a r i a t i o n , and the objective of t h i s chapter i s to es t a b l i s h that the coexistence of species i n such an environment can be dependent on environmental f l u c t u a t i o n . Construction of an elemental model with some claim to b i o l o g i c a l realism i n which environmental f l u c t u a t i o n i s necessary to maintain coexistence i n the absence of any other s t a b i l i z i n g factor w i l l be taken as evidence that in more complete and complex models the temporal pattern of environmental v a r i a t i o n must be considered part of the basis f o r coexistence. U n t i l recently, the p o s s i b i l i t y of persistent multispecies communities maintained by environmental f l u c t u a t i o n had not been investigated through mathematical models. Although Hutchinson[20] suggested the p o s s i b i l i t y of an opportunistic community of phytoplankton which p e r s i s t s when no s i n g l e equilibrium point e x i s t s due to environmental f l u c t u a t i o n . 10 example models were not reported u n t i l Grenny,Bella and Curlf211 and Stewart and Levin[22] described coexistence on a s i n g l e nutrient of two species i n a variable flow chemostat when each species has d i f f e r i n g nonlinear nutrient uptake functions of nutrient concentration. Considering species having non-autonomous nutrient uptake rates, Armstrong and McGhee[23] give a model i n which n species coexist when p a r t i t i o n i n g a s i n g l e conservative resource i n a periodic environment. Their model assumes d i s j o i n t growing seasons of each species with each species dying o f f s u f f i c i e n t l y r a p i d l y a f t e r i t s growing season, releasing the resource so that i t i n t e r f e r e s only s l i g h t l y with the growth of other species. They show that species c h a r a c t e r i s t i c s can be chosen so that a l l species p e r s i s t on a s i n g l e resource i n a homogenous environment. A more comprehensive model w i l l be given i n which a regenerating nonconservative nutrient can support many consumers growing concurrently i n a f l u c t u a t i n g environment, but only one consumer i n a fixed environment. The eguations of the model are conceived as biomass rate of change eguations subject to biomass generation, t r a n s f e r , and l o s s l i m i t s with biomass flow d i r e c t i o n and transfer e f f i c i e n c y c o n s t r a i n t s . The c o e f f i c i e n t s of these equations are assumed functions of a f l u c t u a t i n g environmental variable. The model may be regarded as a generalization to time varying c o e f f i c e n t s of Lotka-Volterra-Terhulst equations representinq a s i n g l e prey with density dependent s p e c i f i c growth rate on which many predators subsist 11 competing only through t h e i r e f f e c t s on t h e i r food supply, the prey. The s i m p l i f i c a t i o n to - a l i n e a r dependence of s p e c i f i c growth rate on species biomass density made i n the L-V-V model i s b i o l o g i c a l l y non-generic, and excludes self-sustained f l u c t u a t i n g solutions of autonomous eguations which have been shown by McGhee and Armstrong!; 2U ] and Koch[ 25 ][ 26 ] to be capable of creating a persistent two predator one prey system. The l i n e a r i t y assumption ensures that f o r f i x e d c o e f f i c i e n t s , only one predator can p e r s i s t on one prey so that the model i s d e t e r m i n i s t i c a l l y unstable. This r e s t r i c t i v e assumption i s made to exclude another p o t e n t i a l basis for coexistence so that the e f f e c t s of environmental f l u c t u a t i o n alone may be studied. By the same l o g i c , the ordinary d i f f e r e n t i a l eguations model i s selected as the most competitive case of a more general p a r t i a l d i f f e r e n t i a l eguation model in which the environment i s s p a t i a l l y heterogenous and each predator moves to those areas where i t s environmental c h a r a c t e r i s t i c s and a v a i l a b l e food supply r e s u l t i n the greatest predator growth. When assumption of homogenous environment, nutrient, and predator densities i s made, mitigation of competition due to s p a t i a l segregation of predators i s precluded along with coexistence a r i s i n g from a v a i l a b i l i t y of s p a t i a l refuges from competition. S i m i l a r l y , assumption of a s i n g l e undifferentiated food supply i n t e n s i f i e s competition over a more r e a l i s t i c model where multiple food supplies and predator preferences among these food supplies can create coexistence. Age structure and other differences within species biomass are ignored, therefore the model w i l l be an exact representation only i f the c h a r a c t e r i s t i c s of d i f f e r e n t i n d i v i d u a l s are equivalent per unit of biomass. This assumption, although inaccurate for species which have d i f f e r e n t l i f e stages, i s more accurate than assuming population numbers as state variables and regui r i n g c h a r a c t e r i s t i c s be equivalent per i n d i v i d u a l . 1.1 Single Nutrient Multiconsumer Varying Environment Model Equations 1.1 and 1.2 together with i n e q u a l i t i e s 1.3 to 1.9 represent the model. Biomass i s qenerated by the nutrient N at a rate rN - bN2 where r i s the i n t r i n s i c rate of nutrient increase and -bN represents i n h i b i t i o n of nutrient s p e c i f i c growth rate by competition for an unspecified resource such as space. Nutrient biomass i s u t i l i z e d by each consumer at a rate Na ix i, where parameter a i # attack rate, can be regarded as the proportion of nutrient which a unit of consumer biomass can u t i l i z e i n a unit of time. The u t i l i z e d nutrient i s assimilated i n t o consumer biomass at an e f f i c i e n c y e^, and consumer biomass i s l o s t at a r e s p i r a t i o n rate ^ The constraint a i(v) > 0 defines x ± as a consumer, the condition m.> 0 l i m i t s the source of x biomass to assimilated 1 i nutrient, and the maximum e f f i c i e n c y e <1 r e f l e c t s losses ' 1 max inherent i n energy conversion. Maximum rate r represents the 3 1 max physical l i m i t a t i o n s on the maximum growth rate of nutrient, a. represents a l i m i t a t i o n on consumer attack rate, minimum 13 i n h i b i t i o n parameter b . ensures that nutrient, even vhen min u n u t i l i z e d cannot grow without bound and i s ultimately limited by some unspecified resource and the condition m. . >0 excludes species which have no re s p i r a t i o n requirement. The remaining constraints are made to simpl i f y analysis by excluding consumers which feed at negative net e f f i c i e n c i e s e i, which have unbounded re s p i r a t i o n rates and nutrients which may not be able to grow even i n the absence of consumers. Hence, we consider the model n dN/dt = N {r(v) - b(v) N - j[ a.(v) x.} i=l 1 1 dx±/dt = x ± {e±(v) a±(v) N - m^v)} i = 1, n 1- 1 1.2 N(0) > 0 1.3 x±(0) > 0 i - 1, n. 1.4 The c o e f f i c i e n t s of equations 1.1 and 1.2 are functions of an environmental variable v{t), a function of time only. The functions have the followinq r e s t r i c t i o n s : rmax ± r ( v ) ± rmin > ° 1 ' 5 b > b(v) > b , > 0 1.6 max — — min a. > a.(v) > 0 imax — i — 1.7 14 1 > e > e. > e. (v) > 0 1-8 max — imax — i ' — m > m. (v) > m. . > m . > 0 1-9 imax — i — imin — min 1.2 General Bounds On Model Biomasses The solutions of the model equations with the given r e s t r i c t i o n s have c e r t a i n general bounds which are proved i n Appendix I. Namely, a l l biomasses have upper bounds, are non-negative f o r po s i t i v e i n i t i a l biomasses, and cannot become zero i n f i n i t e time. Nutrient biomass has an upper bound N = max {r(v)/b(v)> 1.10 independent of i n i t i a l N provided N(0) < N f which i s reasonable since there can be no generation of nutrient biomass by the terms rN -bN 2 i f the bound i s exceeded. T o t a l biomass has an upper bound N(t) + T x. (t) < N + r /4b m 1.12 v ' L i — max min min independent of N {0) and x{0) provided N(0) < N and J" x. (0) < r 2 /4b . m ' „ . 1 - 1 3 v / — L i — max min min This i s the amount of biomass which would be present i f the nutrient were at i t s upper bound and consumers were s t a t i c and being supported by the maximum production of the nutrient at 15 '}.' 10056 e f f i c i e n c y of conversion and minimum r e s p i r a t i o n l o s s . The biomasses must also remain p o s i t i v e for f i n i t e t . That i s , even f o r s t o c h a s t i c v ( t ) # no biomass becomes zero for f i n i t e t , although zero may be approached r a p i d l y . 1.3 Deterministic I n s t a b i l i t y Of Multiconsumer System When the environmental variable v i s a constant, i t can be shown (Appendix II) that at most one consumer may survive i n the system. Consider the nutrient l e v e l necessary to maintain a species at a constant l e v e l by reguiring dx i/dt = 0. From equation 1.2, that maintainance l e v e l i s 1.14 N ± = N±(v) = mi(v)/ei(v) a±(v) The species k for which N =min N i s the only survivor, , since k i i i t can grow more rapidly than any other species, and reduces N to N where every other species has a negative growth rate. k The multiconsumer community i s therefore d e t e r m i n i s t i c a l l y unstable since i t cannot p e r s i s t i n a fixed environment. 1.4 Approximations To A Persistent Multispecies Solution In A-Periodic Environment We w i l l now examine the persistence of a multispecies community by giving error bounds f o r c e r t a i n approximations so that a r e a d i l y calculated approximation to the mean consumer biomasses w i l l be close to the actual s o l u t i o n . 1.4.1 Memoryless Nutrient Approximation I f the dynamics of the nutrient N are rapid enough so that a guasi-equilibrium of nutrient i s c l o s e l y approached before 16 appreciable changes i n the environment or consumer biomasses occur, that equilibrium value may be substituted f o r N in the consumer qrowth equations. Let us examine the errors i n t h i s approximation of the dynamic nutrient biomass by a memoryless expression i n which the nutrient state does not appear. Equation 1.1 may be written as dN/dt = Nb(v) (K(t) - N) 1 - 1 5 where r(v) - I a±(v) x.(t) K(t) = b(v) 1.16 I f c o e f f i c i e n t s b and K were constants, equation 1.15 would be the f a m i l i a r l o g i s t i c growth law i n a rearranged form with so l u t i o n N(t) = - j K^O 1 - 1 7 N(t) = —:— K = 0 1 - 1 8 t(oT + b t For f i x e d b and K, where K > 0, H(t) approaches K uniformly from a l l N(0) > 0. The approach from N(0) > K i s greater than exponential rate bK, and the approach from 0<N(0) <K i s l e s s than exponential rate bK and greater than exponential rate bH(0). I f K<0, N approaches zero more r a p i d l y than equation 17 .1.18.. Defining the average N over i n t e r v a l t-T < x < t as , t <N(t)> - ± J N(x) dx , 1.19 t-T f o r f i x e d b and K>0 ve obtain by integr a t i n g equation 1.17 over t h i s i n t e r v a l * H ^ > - V + 1 fn N ( t ) U 2 0 <N(t)>T - K + — l * ^ j f j . How consider the time varying case with a given K(t) for which 0 < K J < K (t) < K . The region K . < N < K i s an min max ' mxn max at t r a c t o r f o r N since dN/dt < Nb(K - N) 1.21 — max and by integration N(t) < - 1 + <»7n7 " lT—> J, b d T> K T VN(0) K ' max > — 1.22 max max o From equation 1.22, N(t) < K i f 8(0) < K . For N(0) > 3 max max K , Nft) < K as t increases i n d e f i n i t e l y , but dN/dt < 0 max max when N=K and dN/dt= 0 only i f K (t)E K , therefore N(t)<K „ max max max af t e r some s e t t l i n g time. S i m i l a r l y , for K . > 0 min 1.23 dN/dt _> Nb(K m i n - N) and by integration 18 1> N(t) > 7 T Jt ~ 1 7 + (NToy - r r > e x p { - K m i n /„ bd^> 1- 2 « min min 0 thus N (t)>K . a f t e r some s e t t l i n g time. Thus a l l t r a j e c t o r i e s min eventually enter the region 0<K . < N < K__ and any 1 ' mm max ' t r a j e c t o r i e s within the region remain within, hence steady state solutions f o r N are bounded by the region 0<K . < N(t) < K min max The average error i n using K as an approximation to <N{nT)>T when K changes abruptly at f i x e d i n t e r v a l T i s , from eguation 1.20, bounded by 1 N l |<N(nT)>T-K| < max 1^ ^ — | K 2 5 1' 2' Given the bounds on N at steady state and bounds on b (eguation 1.6) we have |<N(nT)>T-K| < r - ± T £ n f -min min as a bound on the average error over an i n t e r v a l T i n length of using K as an approximation for <N(nT)> Consider now a given d i f f e r e n t i a b l e K ( t ) . The actual solution N(t) tracks K(t) i n the sense that i f and only i f N<K i s dB/dt > 0 and i f and only i f N > K i s dN/dt <0, which i s evident by inspection of eguation 1.15. Equality i s attained i f and only i f dN/dt=0. Consider the deviation n(t)=N(t)-K(t) of N from the target K. In a steady state s o l u t i o n , the maximum of a l l l o c a l maxima of n(t) i s the maximum error, and the minimum of a l l l o c a l minima of n(t) i s the minimum er r o r . Both these extrema are 19 included i n the set of a l l possible points i n N where d n/dt=0. Solving f o r those N* for which the error i s stationary, 0 = dn/dt = dN/dt - dK/dt 1-27 and s u b s t i t u t i n g i n Equation 1.15 for dN*/dt gives 0 = N*b(K - N*) - dK/dt Solving t h i s quadratic i n N* y i e l d s 1.28 N* = | ± /(K/2)2 - K/b where K =• dK/dt. 1-29 This set forms two separate r e a l functions of time provided the roots are r e a l , that i s , 4K/bK2 < 1 for a l l t. 1-30 The lower branch remains below any steady state s o l u t i o n N(t) > K , provided min 1.31 K/bK < K for a l l t. min Be now consider only the upper branch, f o r which the values of n*(t) where an extremum can occur are n*(t) = N*(t) - K(t) =. -K/2 + /(K/2)2 - K/b 1 « 3 2 Therefore the maximum and minimum deviations are bounded by n (t) < n* and n(t) > n* . which combine to y i e l d max min 20 |n(t)| < max {K/2(-l + A - 4K/K2b)} 1 - 3 3 which for 4K/K2b « 1 approaches the condition | n ( t ) | < -ainK/Kb. Summarizing fo r a given K(t) d i f f e r e n t i a b l e , i f f o r a l l t 0 < K ' < K(t) < K 1-34 min — — max 4K(t)/b(t)K2(t) < 1 U 3 5 K(t)/b(t)K(t) < 1.36 then the deviation n (t)=N (t)-K(t) of the steady state solution from the target K i s bounded by |n(t)| < max {K/2(-l + A - 4K/K2b)} . 1.37 He have now obtained s u f f i c i e n t conditions to bound the error i n the memoryless nutrient approximation f o r a given continuously varying K (t) (equation 1.37) or the average error f o r a given K(t) changing at discrete times nT (eguation 1.26). 1.4.2 Memoryless Nutrient Consumer System flaking the s u b s t i t u t i o n K(t)=N(t) in equations 1.2 we obtain 1.38 dx±/dt = x ± { e ± ( v ) a ± ( v ) [r(v) - £ a^(v) X j ] / b ( v ) - m ^ C v ) } Defining a l e s s redundant set of functions, l e t 1 39 g±(y) - e ± ( v ) a ± ( v ) r ( v ) / b ( v ) - m ^ v ) 21 u,(v) = a,(v)/b(v) 1.40 w±(v) = e ±(v) a±(v) and logarithmically transforming variables to y i = l n ( x i ) , dx±/dt = x ± {g±(v) - w±(v) I u(v) x } 1.41 1.42 becomes dy_/dt = j> - w u' exp(y_) 1.43 where y.,.g, w,u, exp (v.) are column vectors and * denotes transpose. These are Lotka-Volterra-Verhulst competition eguations with a p a r t i c u l a r l y structured time varying i n t e r a c t i o n matrix composed of the outer product of w<v) representing a s s i m i l a t i o n parameters, and u(v) representing u t i l i z a t i o n parameters of the competition f o r nutrient. The growth vector ij(v) has absorbed both assimilation and r e s p i r a t i o n e f f e c t s . By conditions 1.6,1.7, and 1.8, the elements of the i n t e r a c t i o n matrix C(t) = w(v) u* (v) are non-negative. 1.4.3 Bounds On Consumer Biomass Variation For Periodic Solutions With Periodic Environment Assume a periodic solution x(t) of the same period as the environment. This i s plausible since s u b s t i t u t i o n of an assumed T periodic s o l u t i o n N(t),x{t) of undetermined form i n t o equations 1.1 and 1.2 involves only i n t e g r a l multiples of the fundamental environmental frequency by the following argument. Piecewise d i f f e r e n t i a b l e environmentally dependent c o e f f i c i e n t s denoted by J3(v(t)) may be expressed, for a given piecewise d i f f e r e n t i a b l e v ( t ) , as Fourier s e r i e s i n cot where u i s the 22 fundamental frequency of v (t) . Substitution of Pourier s e r i e s i n ojt with undetermined c o e f f i c i e n t s f o r a T periodic steady state solution of N(t) and x(t) in t o equations 1.1 and 1.2, along with the known series f o r £{v(t)) y i e l d s a set of equations involving harmonics i n cot. Only i n t e g r a l harmonics of a are present since products of the biomasses and parameters i n the r i g h t hand sides of eguations 1.1 and 1.2 r e s u l t in sums whose terms exp(jnojt+jmcjt+jka)t+jiojt) resolve i n t o frequency components which are i n t e g r a l multiples of the fundamental environmental freguency w only. Since a l l frequencies i n the derivative of the undetermined s o l u t i o n must be i n t e g r a l multiples of u, the r e s u l t i n g s o l u t i o n of the harmonic balance equations obtained must be T pe r i o d i c . The absence of subharmonics of to depends on the s p e c i f i c qrowth rates (or per unit biomass qrowth rates) being l i n e a r i n the biomasses. I f , f o r example u t i l i z a t i o n were a saturatinq function of nutrient density, then s u b s t i t u t i o n of an assumed T periodic solution would produce terms which were not i n t e q r a l harmonics of the fundamental, i n v a l i d a t i n q the assumed s o l u t i o n form. This allows periodic solutions not T periodic and non-periodic solutions when s p e c i f i c qrowth rates are not l i n e a r functions of biomass. Each consumer biomass x i has a maximum rate of exponential increase bounded by dx^/dt < x i 0 i where Di=max e i ( v ) a i ( v ) S -• i ( v ) and a minimum rate of exponential decrease bounded by a ^ / d t > - x ± l ± where L ± < max 3 B « i m a i 23 Dropping subscripts, the most variable periodic x(t) f o r a given mean <x> and period T i s obtained when dx/dt=xU f o r part of the period and dx/dt=-xL f o r the remainder of the period, where the time t at which growth switches from maximum increase to maximum decrease, i s chosen t o make x(t) p e r i o d i c . This corresponds, over the period 0 < t < T to the so l u t i o n x(t) = x . exp(Ut) 0 < t < t - LT/(U+L) min s x(t ) = x - x_.._ exp(ULT/U+L)) max mm < t < T X<T> * xmin which has mean given by < x > . x exp(qT) - 1 = T min qT max qT where q = UL/(U + L) In terms of the mean, the extremes of i . are: Ximin = < X i > S*n***iq±T) " « X j_.„ - x ± q^ / C l - «p(-q±T)) imax Therefore the maximum so l u t i o n has a bound 1.44 1.45 1.46 1.47 1 - exp(qT) 1-1*8 1.49 1.50 1.51 deviation from the mean of a periodic 24 q±T - 1 + exp(-qit) |x,(t) - <x > | < <x >d d 1.52 1 i T i i i 1 - exp(-q.t) and x± (t) i s given by 1.53 c ±(t) = <x±> (1 + S ±(t)) |6±<t>| < d ± The r e l a t i v e deviation bound d^, as q « 1 s i m p l i f i e s to approximately d± <q T/2. 1.4.4 Bounds On Biomass Means For A Periodic Memoryless Nutrient Consumer System Assume that v(t) i s periodic with period T and a steady state solution of x{t) r e s u l t s with the same period. Integrating equations 1.43 over one period. 1.54 y(t+T) t+T 0 = / dy = / £(t) - C(t) x(t) dt y(t) t But from equation 1.53, equation 1.54 becomes t+T <£>T - f J c<c) (1 + 0(0) <x> dt U 5 5 Where I i s the i d e n t i t y matrix and D i s a diagonal matrix with diagonal e n t r i e s d = 6^t). Each term of the i n t e g r a l i s t+T n t+T 1 1 r 1 i / c ± j(t) (l+6j(t))<xj>dt o <c±j><xj> + Y <Xj> J t c i j ( t ) 6 j ( t ) d t .56 But since c (t) >0 and |5. (t) | < d. i j J J .t+T -J 57 <cij>dj < ^ / c i j ( t ) 6 j " ( t ) d t - ^ i j " d' J therefore 25 i t + T i / c,.(t) 6.(t) dt = <c,.> e., |e.,| < d. Equation 1.54 then i s qiven by <g> - <C> (I+E) <x> = 0 1 " 5 9 where e.. are the elements of square matrix E. I f <C> i s nonsinqular and defining x*= <C> ~ l <g> as the mean approximation, required to be f e a s i b l e , that i s x1*>0, then from equation 1.59 the mean i s qiven by , 1.60 <x> = (I+E) x* How the inverse may be expanded i n an i n f i n i t e s e r i e s (I+E)'1 = I - E + E 2 - E 3 ... 1 - 6 1 which can be v e r i f i e d by m u l t i p l i c a t i o n . We can therefore qive the mean of the sol u t i o n as <x> - (I - E + E 2 - E 3 ...) x* 1 - 6 2 or <x> = x* + err 1.63 We now wish to bound err . The norm of Ex* s a t i s f i e s the r e l a t i o n ||Ex*||2 = x*'E'E x* < ||x*||2 1.64 where X i s the l a r q e s t eigenvalue of the matrix E * £ r 4 6 1 . max Since the matrix E'E must be p o s i t i v e semidefinite, the eigenvalues are nonneqative. Therefore the trace bounds the aaximum eigenvalue and since T (E'E) 1 n J d 2 1 " 6 5 r i=l 26 ve have that | |Ex*| | <_ p||x* Where P ' y n I d i-1 S i m i l a r l y and hence i s by equation 1.67 and 1.69 bounded by CO l l e r r l l < I 11^*11 £ I M I I P m=l m=l Sunning the geometric s e r i e s ve obtain lle r r l l < I|x*|| p/(l- P) p < 1 1.66 1.67 |E2x*|| < p||Ex*|| < p2||x*|| Emx* < pm||x*|| 1 - 6 9 Therefore 00 err = I (-E)m x* 1.70 m=l m 1.71 1.72 Since l e r r ^ l < Ilerr11 the mean biomass of each consumer i s bounded by <x±> > x ±* - ||err|| . 1 « 7 3 I n t u i t i v e l y t h i s says that ve have no more than enerqy err vhich cannot be accounted for i n the approximate solut i o n x*. 27 I f that amount of energy exceeds the energy l e v e l of the smallest biomass, then coexistence cannot be assured, since the energy may be taken from that biomass producing e x t i n c t i o n . Summarizing, i f the mean i n t e r a c t i o n matrix i t + T 1 t + T 1 7U <C> = ± J C(t)dt = f j w(t) u'(t)dt i s nonsingular and i f the mean approximation x* neglecting consumer v a r i a t i o n given by x* = <c> 1 <g> 1.75 i s f e a s i b l e , that i s x * > 0, then there i s a maximum absolute error i n each term of x* given by eguations 1.72, 1.67 and 1.52 involving the bounds on r e l a t i v e deviations of consumer solutions from t h e i r means. h looser bound on err may be obtained by noting that £ d i 2 < nd 2 vhere i n t h i s case p<nd where d i s the maximum max max max deviation of any species from i t s mean, and n i s the number of species. Heuristic s u f f i c i e n t reguirements f o r persistence of a given stable, f e a s i b l e x* are therefore min x * > ||x*|| nd 7(l-nd ) nd < 1 1.76 i i M ii m a x ' ^ max' max This requires that for solutions with a wide range of mean magnitudes, the r e l a t i v e deviations from the means must be smaller to ensure coexistence than f o r a system with more equitable mean biomasses. 28 1.4.5 Constructive Conditions For A Hodel With & Persistent Solution In order to include the e f f e c t of error r\(t) i n the s u f f i c i e n t conditions for a persistent s o l u t i o n , a set of quadratic rather than l i n e a r i n e q u a l i t i e s must be considered i n order to determine bounds. The bounds on tarqet nutrient K(t) involve the consumer solut i o n bounds and the r e s u l t i n g error bound on n(t) i s not a l i n e a r function of the mean populations,thus the i n e q u a l i t i e s are i n t r a c t a b l e . We s h a l l therefore show that i f a system i s qiven with a periodic v ( t ) , f e a s i b l e and stable mean approximation x*, and non-overexploited nutrient when x=x*, that i s . then a coexistent system can be obtained from t h i s system by s c a l i n q of c o e f f i c i e n t functions., Consider the modified system where constant p o s i t i v e s c a l i n g f a c t o r s s and c multiply the co e f f i c e n t functions as i n K(v,x*) = (r(v) - I a±(v) x±*)/b(v) > 0 for a l l t, 1.77 dN/dt = N(s(r-bN) - £ a±x±) 1.78 dx,/dt = xjLc(eiaiN-mi) 1.79 This system has a mean approximation x* =s x* where x* corresponds to s=1,c=1. The target nutrient at approximate mean consumer l e v e l s , ! 29 1.80 K* « K(v,sx*) = (r(v) - I a±(v) x*)/b(v) i s not a function of s. The maximum nutrient bound N remains unchanged, and raaximum rates of change of x are multiplied by c, therefore the r e l a t i v e deviations d^ccL^ approximately. Therefore the deviation from mean approximation sx* can be made a r b i t r a r i l y small by making c a r b i t r a r i l y small. Then 1 81 K(v,x) K(v,sx*) = K* and the error (t) can be made a r b i t r a r i l y small since |n(t)| < max ^ | {-1+ A - 4K*/sbK*2} 1.82 decreases as s i s increased. Since the two sources of error i n the mean approximation jc* can be made a r b i t r a r i l y small by increasing s and decreasing c, i t w i l l always be possible to construct a persistent system from a f e a s i b l e and stable approximate mean and non-overexploited nutrient system. This i s done by increasing nutrient productivity and decreasing the energy transfer rates of consumers u n t i l a persistent solution i s obtained by numerical solu t i o n of the d i f f e r e n t i a l eguations. 1.4.6 S t a b i l i t y Of A Related Time Invariant System The related time invariant eguations £ = <£> -- <C> exp(y_) 1.83 aay approximate the f u l l dynamics of eguations 1.43 whose steady-state solution has small deviations from a f e a s i b l e x*. The s u f f i c e n t conditions for a f e a s i b l e stationary point x* of 30 in-equations 1.83 to be an asymptotically stable equilibrium i n the la r g e are that the eigenvalues of <C> have a l l r e a l parts positive[ 14 ]. . When the terms be^ are not functions of an environmental variable, the change of variables exp(y i)= be iexp(z i) gives equations = <g> - w exp(z) 1.84 where W = J W(T) W'(T) dx 1.85 Since W must be non-sinqular f o r x* to be unique and / a/ W(T) W' (T) jjtx _> 0 for a l l 1.86 the i n t e r a c t i o n matrix H of equations 1.85 i s p o s i t i v e d e f i n i t e . Therefore under the r e s t r i c t i o n s that be. not vary, a l l I r e l a t e d time invariant equations 1.83 which have a f e a s i b l e stationary point x* are stable, since p o s i t i v e d e f i n i t e matrices have p o s i t i v e r e a l eigenvalues. 1.5 Nutrient Growth, Environmental Fluctuation And Consumer Species Diversity. I f nutrient i s assumed to be i n competition f o r a resource S assumed to be space where each unit of nutrient occupies area a(v) and growth takes place to f i l l the unoccupied space S-a(v)N at the rate i(v), the nutrient growth equation i s dN/dt = N [A(v) (S - o(v)N) - I a ±x ±] 1 - 8 7 implying 31 r(v) = A(v) S and b(v) = A(v) o(v) For equivalent systems d i f f e r i n g only in the magnitude of J!(v) corresponding to d i f f e r i n g s only i n eguation 1-78, the more productive system w i l l be able to support a larger biomass of consumers, and through the more rapid recovery of the nutrient r e f l e c t e d by the decrease i n nutrient tracking error T)(t) , w i l l be able to support a more diverse consumer community. Highly s p e c i a l i z e d consumers growing rapidly i n a small subset only of v(t) w i l l be able to u t i l i z e the nutrient without i n t e r f e r i n g with the subseguent consumers. In the l e s s productive system, e f f e c t s of previous e x p l o i t a t i o n w i l l be adversely f e l t by subseguent consumers, so that highly s p e c i a l i z e d consumer niches cease to e x i s t , r e s u l t i n g i n a l e s s diverse community with broad niche species. The l a t i t u d i n a l gradient of annual s o l a r input i s an i l l u s t r a t i o n of an increase i n £(v) with reference to vegetation growth. The range of environmental v a r i a t i o n such that a multispecies community pe r s i s t s i s i m p l i c i t l y assumed to be such that no one species can be best adapted at a l l environmental st a t e s . This does not necessarily imply a rigourous physical environment, since the physical rigour i s only important i n determining whether a consumer's mean growth rate without competition and at low densities i s p o s i t i v e . The question of 1.88 1.89 32 whether one species can be superior at a l l environmental states i s r e l a t i v e to a l l other a v a i l a b l e competitors,, This i s b a s i c a l l y an empirical question to be answered f o r p a r t i c u l a r cases. He can, however recoqnize that f o r incremental chanqes i n a species* environmental response functions, subject to constraints of current structure, s e l e c t i v e pressure on each species to become unive r s a l l y superior may not e x i s t , even on qeological time scales. Environmental uncertainty, combined with the previously mentioned l o c a l constraints, may make s p e c i a l i z a t i o n i n environmental conditions the predominant d i r e c t i o n of s e l e c t i o n , although t h i s i s not g l o b a l l y optimal f o r a species. I f the resource S i s interpreted as a l i m i t i n g nutrient, such as nitrogen f o r phytoplankton, then an increase in S by enrichment has a q u a l i t a t i v l e y d i f f e r e n t e f f e c t than an increase i n incident l i g h t f l u x i(v). Enrichment a f f e c t s the low density growth rates .g{v) d i f f e r e n t l y f or each species and with consequent disturbances to the x(t) and target nutrient K = (Jl(v)S - I a±x±)/SL(v) a(v) ' 1.90 The r e l a t i v e proportion of species i n the enriched consumer mean approximation may be much d i f f e r e n t from the unperturbed system vheras the increase in l i g h t f l u x merely increases a l l consumers i n the same proportion and leaves standing crop of nutrient unaffected. 33 1,6 Conclusions Assumption of a memoryless nutrient at the time scale of environmental and consumer species change allowed h e u r i s t i c s u f f i c i e n t conditions to be given for a persistent multispecies consumer community subsisting on that nutrient i f a periodic steady state solution i s assumed. Where the nutrient i s not memoryless, constructive s u f f i c i e n t conditions were given f o r a persistent multispecies community. The conditions given are very r e s t r i c t i v e , including the condition that nutrient may never be u t i l i z e d at an i n t e n s i t y r e s u l t i n g i n eventual depletion of the nutrient i f maintained. These r e s t r i c t i o n s were necessary to obtain a so l u t i o n which was mathematically tractable so that i n e q u a l i t i e s involving the solution could be obtained. The tractable solution corresponds to consumer species i n a periodic environment where the numerical response of each biomass over the period of environmental f l u c t u a t i o n i s r e l a t i v e l y small and the changes in nutrient u t i l i z a t i o n are primarily functional responses to environmental change. In the next chapter, we explore numerically the coexistence when environment may be random, numerical response strong, and v e r i f y the assumption of periodic solutions i n a periodic environment. 34 CHAPTER II EXAMPLES OF COEXISTENCE: NUMERICAL SOLUTIONS As a guide to constructing coexistent model communities, natural communities which resembled the given model structure and i t s solutions were sought. A region i s required within which immigration and emigration of biomass i s small with respect to resident biomasses at a l l times. Consumers should be mixed throughout the region without s p a t i a l segregation on the time scale of nutrient biomass variation and should subsist with s i m i l a r discrimination on food sources which compete for the same resource. Regional environmental variables strongly a f f e c t i n g a l l consumer's feeding behavior must be present, and the e f f e c t s of these variables on each consumer's feeding behavior must d i f f e r . The h e u r i s t i c s u f f i c i e n t conditions f o r the approximate so l u t i o n treated i n Chapter I were, too r e s t r i c t i v e for d i r e c t a p p l i c a t i o n , but may i n d i c a t e important features for maintaining coexistence of consumers i n less describable s o l u t i o n s . In the development of the mean approximate s o l u t i o n , i t was reguired that nutrient never be exploited at an un-sustainable rate. Softening t h i s condition somewhat, periods of intense u t i l i z a t i o n of the nutrient, severely depleting the stock! of nutrient followed by a general consumer population crash must be absent. This condition i s met when nutrient i s s u f f i c i e n t l y productive and u t i l i z a t i o n low enough so that density dependent l i m i t a t i o n of nutrient growth i s larqe compared to the net 35 growth rate observed at a l l times. In the examples, only a s c a l a r environmental variable i s nsed although t h i s r e s t r i c t i o n i s not necessary to arguments of Chapter I used to esta b l i s h the approximate s o l u t i o n . This s i m p l i f i e s construction of the examples. Two systems which were selected as guides were S i n c l a i r ' s [27] estimates and data from the Serengeti grassland herbivore community and a simpified c l a s s i c a l lake model of grazing zooplankton feeding on phytoplankton in a temperate lake. Since the major goal of t h i s chapter i s to obtain some example models with persistent solutions, the r e s u l t i n g models are not intended to be accurate quantitative simulations of the natural systems. The most that was attempted was to indicate that coexistent consumer communites with an order of magnitude resemblance to natural systems can be constructed using the given model. The f i r s t example i s a mythical grassland system where consumer rates of u t i l i z a t i o n and r e s p i r a t i o n are slowed to approximately four percent of r e a l i s t i c values. This i d e a l i z a t i o n allows a cl o s e r approach to the s u f f i c i e n t conditions for a v a l i d mean approximate s o l u t i o n . The i d e a l system i s subjected to a periodic yearly environment and the mean approximation compared to the solut i o n . In order to display the robustness of persistence and to test the u t i l i t y of the mean approximation under random environmental conditions, solutions of the i d e a l system were performed f o r a variety of randomly varying environments, displaying q u a l i t a t i v e l y some of 36 the features of persistence. The second example consists of bringing the rates of the i d e a l i z e d example into a r e a l i s t i c range for comparison with the Serengeti herbivore community. The t h i r d example i s a model of phytoplankton-grazer zooplankton making the c l a s s i c a l assumption that phytoplankton form the sole energy source of grazer zooplankton. In t h i s model, coexistence w i l l be demonstrated i n a periodic environment for a system where numerical response of consumers i s large. The s e n s i t i v i t y of the coexistence i s examined under conditions of varying nutrient productivity. 2.1 Example 1: Ideal A s p a t i a l l y homogenous grassland i s assumed where nutrient (grass) maximum rate of increase r i s a l i n e a r function of r a i n f a l l , nutrient density dependence of growth b i s assumed constant, and nutrient q u a l i t y i s assumed proportional to r a i n f a l l . Three consumers u t i l i z e the nutrient. Species 1 i s a lush growth feeder capable of a high attack rate a.±, having a high r e s p i r a t i o n rate m1# but low e f f i c i e n c y of a s s i m i l a t i o n e^ which becomes negative unless g u a l i t y ( r a i n f a l l ) i s high. Species 2 i s a more generalized feeder capable of a medium attack rate a 2 , with a low r e s p i r a t i o n rate, and a medium e f f i c i e n c y of a s s i m i l a t i o n e 2 p o s i t i v e except during low r a i n f a l l . Species 3 has the lowest attack rate a 3 , a low r e s p i r a t i o n rate m^ , and e f f i c i e n c y e 3 high unless r a i n f a l l i s exceeding low. Consumers are assumed to cease feeding, that i s 37 ^=0 i f e f f i c i e n c y i s not p o s i t i v e . Figure 2.1 summarizes the i d e a l model giving the environmental rate functions i n time units of weeks and environmental variable mm/wk of r a i n . This time unit was chosen since i t i s the absolute minimum time scale over which the model i s assumed to be v a l i d i n terms of averaging r a i n f a l l , nutrient q u a l i t y rate change, and mobility of the mythical herbivores over the region. The nutrient parameters are, as w i l l be shown i n example 2, chosen t o approximate the observed productivity and estimated standing crop of green grass i n the long grass region of the Serengeti as a function of r a i n f a l l . The shape of consumer attack c h a r a c t e r i s t i c s were chosen a r b i t r a r i l y to give three d i s t i n c t i v e kinds of consumers as previously described with r a p i d l y calculable c h a r a c t e r i s t i c s . The numerical values were chosen to s a t i s f y three h e u r i s t i c c r i t e r i a for a r a i n f a l l time varying between 8 and 24 mm/wk with a l l i n t e n s i t i e s egually probable within one year. The conditions were f i r s t l y to give a fe a s i b l e mean approximate solution x*, secondly to s a t i s f y a short memory nutrient condition, and t h i r d l y to s a t i s f y a necessary but not s u f f i c i e n t condition for fluctuations from mean biomass l e v e l s to be small i n a random environment. The r a i n f a l l as a function of time i s assumed to have two forms: a periodic t r i a n g u l a r wave of period 52 weeks or a random four season model with each 13 week season's r a i n f a l l assumed constant at a value sampled from a uniform random d i s t r i b u t i o n . a/v) m kg^ wjc aJv) 2^ m kg-wk 0.5 a,(v) 2^ m kg-wk 78.6" 13.3 to mm 2ff v : — o wk e ^ v kg/k, 0.05-^ v:mm/wk 5 kg/kg 0.08 v:mm/wk £ i a 6 v:mm/wk 6 e 5(v) kg/kg| 0.7 l/wk ,0055 i nu(v) l/wk .00551 m 3(v) ) l/wk .0054 FIGURE 2.1 Environmental characteristics of nutrient and consumers for example ideal (to scale). 39 Figure 2.2 i l l u s t r a t e s the two forms which were chosen., The expected freguency (in a s t a t i s t i c a l sense) within an i n t e r v a l Av i s the same f o r both types of environmental s i g n a l , uniform from VJJJJ^ to v m a x and zero outside these l i m i t s . The mean approximation, c a l c u l a t i n g the mean in t e r a c t i o n matrix C and mean growth g as either t+T t+T C = / C(T)dx and g = J g(t)tx or t t 2. 1 v v max max C = / C(v) p(v)dv and g = J g(v) p(v)dv * min min i s the same for both environments. The mean approximation i s x*= 0.031, 0.147, 0.199 which i s f e a s i b l e . The second h e u r i s t i c condition, a short nutrient memory, i s s p e c i f i e d for the random system by requiring that the maximum time constant of the l i n e a r i z e d nutrient eguation about the memoryless nutrient approximation be l e s s than the season length. That i s , assuming that consumers are at the mean approximation l e v e l s x*, the d i f f e r e n t i a l eguation for N i s dN/dt = N(r - bN - \ a±x*) . 2 " 3 Further assuming that N i s very close to K* =(r-^a.±x /b, the memoryless nutrient approximation, the d i f f e r e n t i a l equation f o r N i s dN/dt - bN(K* - N) 2 m i* which l i n e a r i z e d about K* has a time constant(for constant v) of T(v) = l/b(v) K*(v) . 2 m S FIGURE 2.2 Periodic environment of period T and random environment of season T used on example ideal. Both signals have a uniform distribution of probability of r a i n f a l l between vmaz and vmin. 41 The maximum time constant, T =max T (V) was required to be l e s s n v x ' than 13 wk, the season lenqth. For the ide a l system, x =6.75 wk n s a t i s f y i n g the second h e u r i s t i c condition. The bounds on accuracy of the mean approximate solution as derived i n chapter I and therefore the v a l i d i t y of both the above conditions depends on deviations from consumer means being small. For the random environment a necessary but not s u f f i c i e n t condition f o r deviations from the mean being small i s that the deviations over one season be small f o r a l l species over a l l environmental conditions. That i s x (t+T ) t+T Un 1 / , S 1 = 1/ e?(v)a.(v)N-m.dt| < 1 for a l l i , v, 2.6 1 x i(t) 1 1 J t 1 i i and assuming the value of N fixed at K* over the season reguires T C = | l / {e.(v)ai(v)K*(v) - m±(v)}\ > T g . 2 ' 7 The bounds on r e l a t i v e deviation from the means are: d l=.213,d 2 = . 12,d3=. 09 calculated from equations 1.49 and 1.52. For the periodic environment, the c a l c u l a t i o n of error bound on x* for a memoryless nutrient system as qiven i n equations 1.72, 1.67 and 1.52 i s ||err|| = .83 ||x*|| so that the f e a s i b l e mean does not s a t i s f y the s u f f i c i e n t condition f o r the approximation to hold, however these bounds are very loose, so t h i s example w i l l be used despite that. Before proceeding to the f l u c t u a t i n g environments, consider the s i n g l e species e q u i l i b r i a r e s u l t i n g from a fixed v. The maintenance l e v e l s of nutrient f o r each species are N1=.0949 v>18.666,N2=.17 v>13.333,and H3=.2286 v>6.0. Below the indicated thresholds, no l e v e l of nutrient can support the 42 respective biomasses at zero growth. Thus species 1 survives i f v>18.66, species 2 survives i f 18.66>v>13.33 and species 3 survives i f v<13.33 provided the maximum attainable nutrient l e v e l does not f a l l below the maintenance nutrient of the surviving species. In t h i s system, the maximum nutrient l e v e l f a l l s below N at approximately v=11 below which no species survives. The e q u i l i b r i a and time constants of small perturbations from the qiven e q u i l b r i a as obtained by conventional l i n e a r analysis are qiven i n Table I. TABLE I: Fixed Environments and Resulting E q u i l i b r i a Hith Time Constants Fixed Equilibrium Time Constants R a i n f a l l (years) am/wk N x(1) x(2) x(3) x(1) x(2) x (3) 24.0 0.0948 0.234 0.0 0.0 0.33& 6.2 5.3 16.0 0.17 0.0 0.392 0.0 2.0 2.8 9.3 13.3 0.2286 0.0 0.0 0.0074 2.0 2.8 404.0 8-24 var 0.23* 0.031* 0.147* 0.199* * Mean approximation with uniformly varying v 6 Envelope time constant of damped o s c i l l a t o r y response Bote that species 3 has a very narrow range of v i n which i t can dominate, 11<v<13.33 and i s much lower than the mean approximation solution i f v remains fixed at 13.33. Species 1 on the other hand i s much larger at v=24 than the mean approximation and reduces the nutrient to much below i t s maximum sustainable y i e l d l e v e l , .25 and has o s c i l l a t i n g convergence. In the presence of a random disturbance t h i s can produce a sustained f l u c t u a t i n g pest system i n which nutrient i s severely reduced, pest dies back, nutrient recovers, and the process 43 repeats i n a haphazard manner, anable to s e t t l e to the deterministic equilibrium due to weak s t a b i l i z a t i o n ( J e f f r i e s [ 2 8 ] ) . The numerical solutions for f l u c t u a t i n g environments were performed by an error c o n t r o l l e d , v a r i a b l e time step, fourth order Bunge-Kutta integration routine ( OBC RKC ) with a maximum step of 6.5 weeks and r e l a t i v e error tolerance of 0.001. Uniformly d i s t r i b u t e d pseudo-random numbers were generated by the generator UBC FRAND using the same i n i t i a l seed. Note that d i f f e r e n t f i g u r e s are plotted to d i f f e r e n t scales. Figure 2.3 displays the s o l u t i o n to t r i a n g u l a r wave vmax=24, vmin=8 and period 52 weeks. Although the s u f f i c i e n t conditions f o r bounded error i n the mean approximation have not been s a t i s f i e d the mean approximation gives a good order of magnitude approximation to the s o l u t i o n . Convergence to a p e r i o d i c steady state solution i s e s s e n t i a l l y achieved i n 35 years. Figure 2.4 displays the s o l u t i o n to a random environment with the same maxima and minima, and season length 13 weeks. The displayed s o l u t i o n indicates no trend to exclusion of any species over a time of 2000 y e a r s . ; Fluctuations are much larg e r than i n the periodic case. Species 1 now ranges over 10 3 versus l e s s than 10 1, but the mean approximation remains a reasonable guide to the s o l u t i o n . This case was used to t e s t integration accuracy by reducing maximum stepsize by 1/8. No d i s c e r n i b l e difference was evident i n the s o l u t i o n . Figure 2.5 displays the Species 1 (yearly extremes) Species 3 10 -2 ID" 1 H o — t ib" 2. 125yx 0 feyr Species 2 * -J 1 0 " ^ Nutrient (yearly extremes) 10 -2 -1 0 - 1 i o -125yr 0 125yr: FIGURE 2.3 Response of example ideal to triangular r a i n f a l l signal vmax=24, vmin-8, and period T one year. Mean approximation shown as - *. Species 1 10° n i o " 2 -i 10 -4 10° n Species 2 10 -2 10 ,-4 Species 3 10 w -1 10 -2 J 10 -4. 2000yr q 0 10 Nutrient 0 ~~ •'• . " 2000yr° 4 0 ' ' " 2000yr FIGURE 2.4 Response of example ideal to random environment. Season length .25 year, uniform probability distribution between 8-24. Mean approximation shown as 2000yr _ * Species 1 Species 2 10"2H 10 -6 2000yr Species 3 _2 10 ^ H -6 10 Nutrient 0 2000 yr * -i i o " 2 J 10 -6 2000 yr 0 2000 yr FIGURE 2.5 Response of example ideal to random environment with identical conditions to Figur 2.4, but differing random sequence. 47 s o l u t i o n with a d i f f e r e n t random sequence. Recovery from a very large random perturbation at around t=1200yr demonstates a c h a r a c t e r i s t i c time of response of the order of 100yr for the s o l u t i o n . In order to test the r o l e of the memoryless nutrient condition, the season length was shortened to 1.75 weeks (much l e s s than the maximum nutrient memory of 6.75 weeks). For t h i s example int e g r a t i o n maximum step was 1.75 weeks. Figure 2.6 displays the r e s u l t . Species 1 and 2 show a c l e a r trend to e x t i n c t i o n . Once a species becomes i n s i g n i f i c a n t with respect to nutrient u t i l i z a t i o n , i t s growth i s l i m i t e d t o t a l l y by other species and i f i t i s de c l i n i n g , no further decline can help i t s s p e c i f i c growth rate. Species 3 contributes a3X3=.144 u t i l i z a t i o n impact over the whole v range, species 2 contributes a maximum 1.x10—6 u t i l i z a t i o n impact and species 1 contributes 3x10-* u t i l i z a t i o n impact. Since further declines i n species 1 and 2 cannot s i g n i f i c a n t l y increase food supply by reduction of u t i l i z a t i o n , they s h a l l continue to. decline at the evident exponential rates., A g u a l i t a t i v e explanation of the success of species 3 i s that i t i s able to u t i l i z e the standing crop of nutrient without competition during low r a i n f a l l but does not depend on production of the nutrient during these periods very much, since large crops are inherited from previous wetter seasons. Species 1 and 2 i n h e r i t during wetter seasons a nutrient l e v e l which does not have s u f f i c i e n t time to grow up to the high l e v e l s they can s u c c e s s f u l l y e x p l o i t . Note the Species 1 Species 3 FIGURE 2.6 Response of example ideal with fast environment. Season length 1.75 weeks, normal climate 8-24 mm/wk. Exclusion i s rapid (compare time scales to Figure 2.5). 49 generally lower nutrient peaks i n Figure 2.6 as compared to Figure 2.5 i n support of t h i s explanation. The s p e c i a l i z e d species have been excluded due to a long nutrient memory. The e f f e c t of lengthening seasons to 52 weeks i s displayed i n Figure 2.7. As expected the l e v e l of va r i a t i o n of populations has increased. Species 1 has a range of 10 6 as compared to 10 3 f o r the 13 week season. The populations do not show any clear trends to extincti o n but species 3 i s showing increasing variance which may ind i c a t e a solution with ever-increasing variance, implying a decreasing mean. Since the l e v e l s of var i a t i o n i n species are approaching those which would re s u l t in extinct i o n due to minimum v i a b l e biomass anyway, the s t a b i l i t y of the mean was not further investigated. Returning to the 13 week period, we now consider the e f f e c t s of a s h i f t in climate. I f the o r i g i n a l uniform d i s t r i b u t i o n from 8-24 mm/wk i s considered NORHAL climate, we term the d i s t r i b u t i o n 6-24 DRYER and 8-28 BETTER. Since s h i f t s i n the d i s t r i b u t i o n a l t e r the average growth and i n t e r a c t i o n , the mean approximation and conditions must be recalculated. For the DRYER climate a f e a s i b l e mean approximation of x*=0.029,0.0527,0.2469 r e s u l t s . The memory time of the nutrient i s T =11.7 wk and the minimum response time of consumers i s 36 n wk. For a season length of 13 wk, the maximum memory time approaches the season length so that the nutrient may not be considered memoryless. The solution for a season of 13 i s given i n Figure 2.8. Species 2 shows a c l e a r trend to e x t i n c t i o n . apecies i 1 0 % Species 2 10 - 4 J 10 -8 0 Species 3 0' 10 10 -8 2000yr 0 10 2000yr Nutrient 10 -4. 10 -8 20 iJoyr 0 T FIGURE 2.7 Response of example ideal to normal climate with one year season^ Mean approximation shown as -*. *e*&on ±g, VJI O 0 Species 1 10 -i 10 -5. .-10. 10 - i Species 2 10" Species 3 200'0yr 0 I 10 _ Nutrient 10 -5 -10 2000yr 6 2060yr FIGURE 2.8 Response of system ideal to DRYER climate. Season length .25 year, 52 Although i t may appear paradoxical that species 1, favored by wet weather survives while species 2 succumbs i n a dryer climate, t h i s i s e x p l i c a b l e by the consideration that species 2 and species 3 are i n competition much more strongly than species 1 and species 3 which can be v e r i f i e d by examination of the inverse mean i n t e r a c t i o n matrix,, Species 2 i s outcompeted i n i t s most e f f e c t i v e range 13.33-18.66 by an increased species 3 population. Aside from the exclusion of species 2 which i s more a product of the nutrient memory condition, the reduction i n mean approximation of species 2 i s much larger than the reduction i n species 1 even i f coexistence i s obtained. In order to s a t i s f y the memoryless nutrient condition, a season of 26 was used. Figure 2.9 shows the r e s u l t i n g s o l u t i o n . Despite a very large deviation of x^(t) , there i s no apparent trend to exclusion or increasing variance. At the scale shown, the mean approximation appears to have some accuaracy, but the magnitude of deviations makes visu a l estimation of mean d i f f i c u l t on a logarithmic scale. For a wetter climate a f e a s i b l e mean approximation of x*=0.139,0.117,0.011 i s calculated with x =7.6 x =41- The n c s o l u t i o n i n Figure 2.10 displays apparent but very slow exclusion of species 3. The f a i l u r e of the h e u r i s t i c conditions i s a t tributed to the errors i n using x* to c a l c u l a t e T and T . n c As can be seen, x^(t) exceeds x^* s i g n i f i c a n t l y and i s the dominant u t i l i z e r . In addition, i t has a very high attack rate. From the graph, the peak u t i l i z a t i o n impact i s .5 f o r v over 1 Q 0 Species 1 10 -10 J 10 10 - 2 0 1 Q 0 ^ Species 2 * -Ml ID" 1 0 H 10 - 2 0 0. Species 3 10 -10J i p - 2 Q 2o60yr 0 2000yr , ~0 Nutrient 10 — i l O " 1 0 ^ - ^ f 2 0 - f 2000yr 0 20o6yr FIGURE 2.9 Response of system ideal to DRYER climate, season length .5 year. Mean approximation shown as *. Species 1 10 0 _ io-' Species 2_° 10° " 10 -5 4000yr Species 3 Nutrient 10° 4000yr 10 -5 0 0 4000yr FIGURE 2.10 Response of system ideal to WETTER climate, season length .25 yr. Mean approximation shown as *. 4000yr 55 18.66. But unless v i s greater than 20, at which point r=.5, the nutrient w i l l be overloaded, that i s , exploited at above i t s capacity and w i l l be driven to low values. Sharp declines i n nutrient to r e l a t i v e l y low l e v e l s are evident in Figure 2.10 i n comparison to NORMAL (Figures 2.4 and 2.5). These declines and subsequent slow recoveries over a number of seasons would be s u f f i c i e n t to deny species 3 the nutrient l e v e l s assumed in the memoryless nutrient approximation and de l i v e r to the already low predicted biomass the coup de qrace. The foregoing solutions have established the existence of persistent solutions to the proposed model. The u t i l i t y of the mean approximation was v e r i f i e d i n i t s use to construct examples. Persistence was shown to be robust i n the sense that both a periodic and random environment of l i m i t e d bandwidth and amplitude can create persistence in the model. The e f f e c t s of a long memory nutrient, overloaded nutrient, and time scale of environmental v a r i a t i o n on the accuracy of the approximate soluti o n were tested. The h e u r i s t i c conditions were useful i n a n t i c i p a t i n g and explaining some departures from the predicted mean solution, but are not on the whole r e l i a b l e for accurately predicting the persistence of solutions. , 2.2 Example 2: Grassland Bates of example 1 w i l l now be scaled up to make rates of exchange between nutrient and consumers and consumer m o r t a l i t i e s commensurate with the data and estimates of S i n c l a i r ^ 2 7 ] for the Serengeti p l a i n s . The attack rates and mortality rates of 56 example one were multipied by 30 times and the e f f i c i e n c y of species 2 s l i g h t l y increased to .085 to obtain coexistence i n a 52 week period t r i a n g u l a r r a i n f a l l with minimum 8mm/wk and maximum 24mm/wk. The mean approximation recalculated with the changed functions i s x* = 0.00044,0.0097,0.00235 kg/m2, and x =5.75, x =1.35. n c Attack rates were calcul a t e d from S i n c l a i r ' s data f27] for the long grassland consumers at highest consumer densities assuming dry weight of mammals to be .1 of body weight and grass standing crop of .25 kg/m2. The biomasses and attack rates are: ungulates, biomass 0.0057 kg/m2 and attack rate 1.0 m2/kg-day, small mammals, biomass 0.000187 kg/m2 and attack rate 6.1.. m2/kg-day, and grasshoppers, biomass 0.0012 kg/m2 with attack rate 6.3 m2/kg-day. Figure 2.11 gives model c h a r a c t e r i s t i c s . Model rates have been changed to days as a more standard e c o l o g i c a l time unit. Figure 2.12 shows S i n c l a i r ' s grass production data obtained by c l i p p i n g selected grass exclosure plots to obtain dry weiqht of grass production compared to one quarter the monthly maximum sustained y i e l d l e v e l of model grass production with S i n c l a i r ' s mean r a i n f a l l figures per month. Figure 2.13 compares model steady state standing crop, removing 450 kg/ha-mo constant u t i l i z a t i o n , with S i n c l a i r ' s estimated standing crop of green grass i n the long grassland which he obtained by accumulating production estimates and subtracting estimated u t i l i z a t i o n . Since standing crop i s not r e l a t e d to productivity of grass i n a 2(v) 2 m /kg-day m /kg-day e 2(v) kg/kg! 18.6 m^v) l/dayf 0.041 l/day| 2.141 0.085J 0.0291 a-(v) 2 ^ m /kg-dayi 6 13.3 1.2-1 0 °6 e_(v) kg^k, 0.1 0 13.3 m3(v) l/day| 0.027-0' FIGURE 2.11 Environmental characteristics of nutrient and consumers for example grassland, (to scale). GRASS 2000 -r MONTHLY PRODUCTION kg/ha 10004 Long grassland mean production - — © — One quarter model maximum sustainable production I i 1 i 1 1 r- r-N D J F M A M J tr*-J A ~0 MONTH FIGURE 2.12 Comparison of production data (Sinclair) with model production under same mean r a i n f a l l pattern. GRASS STANDING CROP kg/he 5000—i Estimates Model - - © — f i - i 1 1 — — r N D J F M A M J J A S 0 MONTH FIGURE 2.13 Comparison of estimated standing crop of grass on long grassland with model steady state standing crop under 450 kg/ha-mo constant u t i l i z a t i o n with measured mean r a i n f a l l . 59 S i n c l a i r * s treatment, there i s some l a t i t u d e i n choice of l o g i s t i c parameters. For the assumed model, the maximum growth ra t e per day of grass at 24 mm/wk r a i n f a l l i s r=0.0857 per day and the minimum maximum growth rate i s 0.02 86 at 8mm/wk. The assumption of variable quality of grass i s consistent with S i n c l a i r and Braun*s[29] measurements of grass protein content as a function of r a i n f a l l . At a low protein content, grass ceases to supply ungulate nitrogen requirements and growth s u f f e r s , but feeding continues t o supply energy requirements. These measurements were not made on small mammals and grasshoppers, but a l l herbivores preferred green grass, which could indicate a s i m i l a r s i t u a t i o n . The assumption i n the model of increased u t i l i z a t i o n at increased nutrient density disagrees with the observation that u t i l i z a t i o n tends to be constant but qua l i t y v a r i a b l e . The model makes no correction for quality except i m p l i c i t l y i n attack functions which switch o f f at some qu a l i t y . However the e f f e c t s on consumers are s i m i l a r ; a higher r a i n f a l l produces an increase i n consumer growth. The coexistent solution f o r a triangular r a i n f a l l period of 52 weeks varying between 24 mm/wk and 8mm/wk i s given i n Figure 2.14. The mean approximate s o l u t i o n i s some guide to the solu t i o n and l e v e l s of the observed ungulate, small mammal and grasshoppers are comparable i n magnitude to species 3, species 2, and species 1 res p e c t i v e l y . The triangu l a r wave i s somewhat more equable a seasonal pattern than the observed mean monthly extremes of 3.7mm/wk (July) and 35 mm/wk (April) . Species 1 (yearly extremes) kg/m 10 10 -7 0 Species 2 (yearly extemes) kg/m2 K f 3 J SM -A Species., 3 kg/m 10 M 10 ~60yr 6 it lb: yr Nutrient 2 kg/m J flfffifiltw 10~3J 60yr Jyr FIGURE 2.14 Response of system grassland to triangular rainfall signal, period one year. Starting time is at maximum rainfall of 24 mm/wk. 61 Coexistence has been demonstrated in a model of comparable time, production, and biomass l e v e l s to a natural system. This i n no way validates the example model as a representation of the natural system since dynamic data were i n s u f f i c i e n t to test d e tailed assumptions of the model, p a r t i c u l a r l y environmental response function shapes and l e v e l s . The higher l e v e l s of v a r i a b i l i t y i n the model would i n a more complete model be lessened by segregation of consumers i n space, s p e c i a l i z e d adaptations to food shortage such as resting l i f e stages, energy storage and migration and saturation of u t i l i z a t i o n and a s s i m i l a t i o n under conditions of food abundance. Inclusion of these factors would not have allowed the role of a f l u c t u a t i n g environment and varying responses to environment to be c l e a r l y demonstrated as a basis f o r stable coexistence on a single resource. 2.3 Example 3: Lake This example i s an elemental model of pelagic plankton i n an o l i g o t r o p h i c lake assuming food generated within the pelagic region. Nutrient f o r f i l t e r feeding zooplankton i s i d e n t i f i e d with a phytoplankton food source whose primary productivity i s a l i n e a r function of average lake temperature. This f u n c t i o n a l dependence of productivity i s consistent with observed l i n e a r dependence of saturation l i g h t i n t e n s i t i e s of generalized photosynthesis rates up to 30 deg.C. [30 J . Since l i g h t i n t e n s i t y at lake surface over daylight hours i s t y p i c a l l y above saturation photosynthesis l e v e l s [31j the assumption of l i n e a r 62 temperature dependence of productivity i s p l a u s i b l e . The f i l t e r feeding zooplankters of temperate lakes have f i l t e r i n g rates which are strongly dependent on water temperature and which are roughly constant per unit of biomass within species (Wetzel [32] Buckingham[ 33 ]) . Saturation of feeding rate at higher food density i s neglected i n t h i s model for the previously argued purpose of s i m p l i c i t y , and i t may be argued that experimental evidence of saturation at time scales of minutes i n laboratory studies i s not relevant i n the longer time scales of days and weeks due to a s s i m i l a t i o n adaptation [34]. Ranges f o r standing crop and production of phytoplankton [35] i n o l i g o t r o p h i c lakes imply ranges of r/b= .2 - 1.0 g/m2 and r= 0.01 - 3.0 day. Ranges for zooplankton f i l t e r i n g rates[36], e f f i c i e n c y of assimilationf37] and biomass turnover rates [38] imply ranges of a=0.05 - .4 m2/g-day, e=.1 - .3, and m= .04 - .1 1/day.. Figure 2.15 gives the example functions. The environmental v a r i a t i o n was assumed sinusoidal with period 365 days. The example assumed two f i l t e r feeders with d i f f e r i n g attack rates as a function of average lake temperature. Species COLD has a guadratic attack rate with maximum at mid range temperature and zero attack rate at high, and low extremes of temperature v a r i a t i o n . Species WARH has a symmetric monotone cubic attack rate with a maximum at maximum temperature and minimum at minimum temperature. The purpose of t h i s example was to i n v e s t i g a t e the 63 64 properties of a coexistent community when the mean approximation solu t i o n cannot be used due to large v a r i a t i o n s i n a l l consumers. The species i n t h i s example, r e l a t i v e to the previous examples were chosen with l e s s overlapping attack functions so that the coexistence would be l e s s s e n s i s t i v e to consumer c h a r a c t e r i s t i c s . The values of e were selected to be .1 . The r e s p i r a t i o n m was taken as varying l i n e a r l y from 0.05 at minimum temperature to .1 at maximum temperature. The i n i t i a l amax=5.0 for both species was chosen to give a sin g l e species equilibrium with comparable phytoplankton and zooplankton l e v e l s . Computations were performed by error c o n t r o l l e d Runge-Kutta integration with a maximum stepsize of 1/1000 of the displayed time and maximum r e l a t i v e error 0.001. Figure 2.16 displays the sol u t i o n with rmax=1.0 b=1.0 a c =a^ =5.0 and the previously qiven values of m and e max max i d e n t i c a l for both species. A sta b l e periodic solution i s e s s e n t i a l l y reached i n one year. Reducing the attack rate of WARM to 2.5 maximum s t i l l produced coexistence as shown in Figure 2.17. Species WARM has a shorter summer peak but higher value, and species COLD has a more pronounced f a l l peak. To t a l zooplankton would show a spring and f a l l peak. Reducing the productivity of t h i s system by lowering r to •1 and b to .1 produced the s o l u t i o n shown i n Figure 2.18. Peak populations have been lowered by approximately .1 and nutrient l e v e l s remain s i m i l a r but more smooth. Further lowering of r Nutrient FIGURE 2.16s Coexistence in example Lake, rmax=1.0, b=1.0, a_ =a„ - 5.0 . Cmax T/max Nutrient Species WARM 10 rl4 10 -7 10 •-3J l O " 7 ^ 0 7 yr Species COLD FIGURE 2.17 : Coexistence in example Lake. Species WARM is a less effective utilizer of nutrient, rmax= 1.0, b= 1.0, aQ m a x = 5.0, max = 2.5 10 .-3 J 10 -7 Ch Nutrient 10 « Species WARM FIGURE 2.18 « Coexistence in a less productive Lake; max rmax= 0.1, b= 0.1, a~max= 5.0, = 2.5 . 7 yr 68 and b to 0.01 caused exclusion of species COLD. The e f f e c t s of reducing r and b are the same as the e f f e c t on a mean approximation solu t i o n as predicted by eguation 1.78 f f . The general l e v e l of nutrient i s unchanged and population l e v e l s are lowered by the same amount. This example i l l u s t r a t e d rapid s t a b i l i z a t i o n to a periodic steady state s o l u t i o n i n a periodic environment. I t also i l l u s t r a t e s the ro l e of nutrient productivity as an element i n the memoryless nutrient condition permitting coexistence i n a flu c t u a t i n g environment. The scales of var i a t i o n and magnitudes i n t h i s example are roughly comparable to observed data, but det a i l e d comparisons w i l l not be made i n t h i s chapter. 2.4 Conclusions The prediction of stable periodic s o l u t i o n s with the same period as the environment was v e r i f i e d f or the model structure under a periodic environment for cases approaching the s u f f i c i e n t conditions of Chapter 1 and f o r examples where biomasses have large v a r i a t i o n s . Coexistence was produced in an example under both periodic and random environments i l l u s t r a t i n g the robustness of persistence of the environmentally s t a b i l i z e d system. The u t i l i t y of the mean approximate s o l u t i o n was i l l u s t r a t e d even when the h e u r i s t i c s u f f i c i e n t conditions were not s a t i s f i e d . The r o l e of a short memory nutrient was demonstrated to be ce n t r a l i n maintaining the persistence of a multispecies 69 community and i n determining the permissible d i v e r s i t y of consumers i n that community. The proposed model was found to be capable of producing coexistence on a sing l e nutrient of multiple consumers i n a s c a l a r environment with magnitude and time scales of the order of two very d i f f e r e n t natural systems. These r e s u l t s suggested that a more serious attempt to model a natural s i t u a t i o n should be made so that the t h e o r e t i c a l abstractions are brought into more intimate contact with actual f i e l d data. Construction of a simulation model incorporating detailed s p a t i a l , environmental and feeding responses of animals was eschewed for two reasons. F i r s t l y , i t was f e l t that such an undertaking was beyond the scope of t h i s t h e s i s and the knowledge of the author. Furthermore i t was f e l t that much of the detailed behavioral and environmental response data would be unavailable so that such a construction would be an exercise i n imagination. The second objection i s that such an approach would not make use of the comparative s i m p l i c i t y of the model unless the model were i n f e r r e d from the the more complete simulation model dynamics. This could be done either deductively, by reducing the simulation by approximation and i d e a l i z a t i o n s or empirically by f i t t i n g the simple model's input and output c h a r a c t e r i s t i c s to the simulation model's input-output c h a r a c t e r i s t i c s . I f the l a t t e r approach i s f e a s i b l e , the simulation model i s redundant since one might as well attempt to 70 f i t the model c h a r a c t e r i s t i c s d i r e c t l y to f i e l d data. The following chapters w i l l describe the methods and r e s u l t s of f i t t i n g model c h a r a c t e r i s t i c s to f i e l d data from a small temperate lake for the four dominant herbivores. 71 CHAPTER I I I HODEL FITTING PROCEDURE The previous chapters emphasized persistence of a model community of known structure, environmental functions of consumers and nutrient, and environmental variables. The r e l a t i o n between the given conditions and persistence i n the model i s i n general g u a l i t a t i v e . Except f o r i d e a l i z e d cases, numerical solution of the model eguations would appear necessary to determine the persistence properties of such a community. Prediction of persistence of a given set of consumers qiven only information about t h e i r presence or absence in habitats with given environmental f l u c t u a t i o n patterns would appear to reguire a general and a n a l y t i c r e l a t i o n between model structure, environment and persistence together with an adequate sample of a l l possible environmental conditions and r e s u l t i n g species incidences. Since neither of these conditions was thought achievable, the general guestion of whether and to what extent environmental f l u c t u a t i o n s t a b i l i z e s a community w i l l be replaced with a more p a r t i c u l a r and r e s t r i c t e d guestion. Can the given model structure comprehend the time dynamics of a persistent natural community? The property of persistence of the model i s i m p l i c i t , i t s sole function being to e s t a b l i s h that persistence of the given model i s possible and therefore i t i s a possible model f o r a persistent natural community., Since a li m i t e d range of possible environmental f l u c t u a t i o n s and concurrent biomass observations w i l l be obtained i n f i e l d 72 measurements, f u l l matching of persistence properties w i l l not be v e r i f i a b l e . The v a l i d a t i o n of a f i t t e d model must therefore r e s t upon the accuracy of independent predictions of biomass l e v e l s from a l i m i t e d sample of environmental observations and conformity of estimated environmental response functions with independently measured components of these functions. The method used to determine the set of environmental functions r (v) , b (v) , (a (v) ,e i (v) ,m ±(v) , 1= 1, n) i s somewhat more than a curve f i t t i n g procedure and l e s s than a s t a t i s t i c a l l y v a l i d parameter estimation. Model d i f f e r e n t i a l eguations are s a t i s f i e d over the time range of given observations and environmental functions are chosen to minimize the deviations of model biomasses from observed biomasses. Prediction of biomass l e v e l s , p a r t i c u l a r l y for a r e s t r i c t e d c l a s s of seasonal flu c t u a t i o n s , may be possible without a uniguely determined functions estimate but for meaningful comparison with independently obtained environmental functions a unique description i s necessary. I t w i l l be shown that the following conditions are s u f f i c i e n t that the system be i d e n t i f i a b l e i n the sense that functions be c o r r e c t l y determined from observations. 1. The system to be i d e n t i f i e d a c t u a l l y has the model structure assumed. 2. The environmental s i q n a l p e r s i s t e n t l y excites the system so that a l l areas of st a t e and environmental space are sampled. 3. Biomasses and environmental s i q n a l are sampled 73 s u f f i c i e n t l y rapidly so that the time variations are characterized. < 4. There must be no measurement errors, that i s observations are exact. Some of these conditions have been rigourously posed and simultaneous state and parameter estimation r e s u l t s obtained for time invariant l i n e a r d i f f e r e n t i a l and difference equations when state observations are corrupted by additive gaussian noise under r e s t r i c t i v e assumptions, but these r e s u l t s are not d i r e c t l y applicable to t h i s model, nor i s a general method av a i l a b l e even for these cases. The guestion of i n f e r r i n g model structure and noise properties from noisy observations remains open even for these i d e a l i z e d systems. The data which were to be used to f i t the model comprise observations of average lake temperature, phytoplankton biomass densities sampled at an in t e r v a l s of approximately one week, and four grazing zooplankter biomasses sampled at approximately two week in t e r v a l s from May 2 to December 12, 1974 i n Eunice lake, a small oligotrophic temperate lake. These data w i l l be more f u l l y described i n the next chapter, but none of the previous conditions are actually s a t i s f i e d by them. The b i o l o g i c a l system almost c e r t a i n l y does not exactly f i t the model structure due to s p a t i a l v a r i a t i o n s in phytoplankton density and consumers, d i f f e r e n t l i f e stages of consumers, and gross s i m p l i f i c a t i o n s of phytoplankton growth c h a r a c t e r i s t i c s . There i s some doubt that phytoplankton form the sole food source 74 of the grazers as well, and predation on consumers not a function of consumer density and lake temperature i s neqlected. The lake temperature has a predominant one year periodic component correlated with biomass v a r i a t i o n . Sampling rate, p a r t i c u l a r l y of the nutrient i s probably i n s u f f i c i e n t to capture rapid changes i n nutrient, since fourfold changes i n phytoplankton densities have been observed i n a one week i n t e r v a l in t h i s lake. F i n a l l y , estimates o f biomass have e r r o r of at least n from In (biomass o b s e r v e d ) =ln (biomass a c t u a l ) + n where n i s normally d i s t r i b u t e d with a standard deviation of .35 since t h i s was the estimated v a r i a t i o n i n weight-length r e l a t i o n s used to cal c u l a t e biomass data. Lacking a developed mathematical approach to these discrepancies, the following approach was adopted. The f i r s t system of Example 3, Chapter II was taken as a t e s t system of known structure on which the f i t t i n g procedure's performance was tested. The e f f e c t s of sampling rates and plausible uncorrelated observation noise levelsi comparable to the f i e l d data on the f i t t i n g procedure were explored. In t h i s manner, the i d e n t i f i a b i l i t y of the model structure i s v e r i f i e d and problems which might be expected i n the planned a p p l i c a t i o n were antic i p a t e d . 3 .t I d e n t i f i c a t i o n Method 3.1.1 Environmental Function Modelling Each unknown function was represented as an orthogonal 75 expansion i n Tchebycheff polynomials with unknown c o e f f i c i e n t s . This representation allows any d i f f e r e n t i a b l e function to be exactly represented over a f i n i t e range of the argument v since the Tchebycheff polynomials form a complete set of orthonormal functions over the range -1<v<1. The Tchebycheff polynomials have the property of maximum convergence [39] inasmuch as they reguire the smallest number of terms of any u l t r a s p h e r i c a l polynomial expansion to achieve an approximation to some given f(v) which does not deviate from the given f(v) by more than a given ± e a t any point within the range This choice of polynomial approximation has a number of advantages computationally. The Tchebycheff polynomials never exceed 1 i n magnitude within the r e s t r i c t e d range of v, therefore c o e f f i c i e n t s used to represent a function are of comparable magnitude, si m p l i f y i n g conditioning problems. Orthogonality of the polynomials eliminates i n t e r a c t i o n of c o e f f i c i e n t s used to determine the functions which reduces i n t e r a c t i o n of parameters and therefore the complexity of the minimization problem within which the functions are embedded. F i n a l l y , the Tchebycheff polynomials are simply calculated by the d e f i n i t i o n s T =1 T=v and the recursion T {v)= 2vT (v) -o 1 n+1 n Tn-l< V>-As an example, the function r(v) i s represented as r(v) 7 r T (v) where R i s nFhe maximum number of terms to be used i n the approximation. The function r(v) i s therefore s p e c i f i e d by R*1 c o e f f i c i e n t s of i t s Tchebycheff expansion. These and • • i 76 c o e f f i c i e n t s of the other functions s h a l l henceforth be referred to c o l l e c t i v e l y as the parameters of the system. In p r i n c i p l e any a n a l y t i c function of r e s t r i c t e d range can be represented but in practice the number of terms in the expansion of the function must be r e s t r i c t e d by the a v a i l a b l e data. This r e s t r i c t i o n aside, t h i s choice of function modelling presupposes nothing about the unknown function. 3.1.2 Continuous Environment Modelling Since the model i s continuous i n time, a means of i n t e r p o l a t i n g between environmental measurements i s necessary to define intermediate values of v ( t ) . A cubic spli n e [40] was f i t t e d through environmental observations using the program OBC SPLINE with observations prescaled within a (-1,*1) i n t e r v a l . The cubic s p l i n e defines a twice d i f f e r e n t i a b l e function of time having minimum curvature passing through a l l observations. 3.1.3 Minimization Problem The solution of the biomass d i f f e r e n t i a l equations for given i n i t i a l conditions x (t^) , given environment v ( t ) , and a p a r t i c u l a r choice of parameters p i s denoted x(t,v,p). Observations of biomass x i=0,n where n i s the number of i k consumers and i=0 corresponds to nutrient are given at times t ^ , k=1,K where the incidence of observations I (i,k) has the value 1 i f species i i s observed at time t, and zero otherwise. An k 'observation* of zero biomass i s considered to be not observed. The t o t a l number of biomass observations i s therefore 77 K n N = I I I(i,k) 3.2 k=l i=0 The objective function of parameters F(P) = | I I (Anx ±(t k) - *nx ± k) 2 k=l i=0 3.3 i s the mean square deviation of the logarithm of the r a t i o of minimum of zero when a l l observed biomasses are i d e n t i c a l to model biomasses. The quantity exp( /F (p)) i s the root mean square r a t i o error frmsre) and gives a r e a d i l y interpreted f i g u r e of merit f o r the error in f i t . For example, i f F{p)=.1, then the rmsre i s 1.37 which i s the error which would be obtained i f a l l model biomasses exceeded observations by 37% or f e l l short of observations by 1/1.37=27%. The i d e n t i f i c a t i o n or estimation problem i s therefore to determine a set of parameters for which F(p) i s a minimum. 3.1.4 Solution Of The Minimization Problem The minimization of F(p) was performed by the Quasi—Newton Variable Metric unconstrained minimization algorithm UBC FMIN written by E. Fletcher and associates of Atomic Energy Research Establishment Harwell, England[41], This implementation of Fletcher's algorithmf42] stores a p o s i t i v e d e f i n i t e approximation to the Hessian ( matrix whose elements are h^= d 2F (p)/dp^p^ ) of the objective function i n LDL* factored form[43] avoiding problems of unbounded matrix e n t r i e s when converging to the minimum of a function having a singular model biomass over observed biomass. This function has a 78 Hessian encountered when the inverse Hessian i s approximated as i n the o r i g i n a l Davidon-Fletcher-Powell method[44"). The minimization procedure used was a general method and did not make use of the l e a s t sguares structure of the chosen objective function. This allows both for improvements i n computational e f f i c i e n c y through e x p l o i t a t i o n of the l e a s t -squares structure or use of objective functions of more complex structure to better achieve the purpose of the modellinq exercise. The l i m i t of improvement on l e a s t sguares objective functions, achievable when state solutions are l i n e a r functions of parameters i s the difference in convergence speed of a Newton-Raphson versus Quasi-Newton algorithm, a factor of n where n i s the number of unknown parameters. This improvement i s not guaranteed even for least sguares problems, and examples have been given[51] where l e a s t sguares algorithms are l e s s e f f i c i e n t than a general algorithm on l e a s t sguares problems. This was a factor i n the choice of minimization algorithm. Gradients of F(p) reguired by Quasi-Newton methods are calculated a n a l y t i c a l l y from the solu t i o n of the state s e n s i t i v i t y equations obtained by d i f f e r e n t i a t i n g the state equations. Let equations 1.1 and 1.2 be denoted as 3.4 D i f f e r e n t i a t i n q with respect t o parameter Pj, the s e n s i t i v i t y d i f f e r e n t i a l equations are qiven by A d z-r J dz. 3 g n 3 g d z d z 3.5 d_ /_!) „ _d_ = + I _ A and - 7 - ^ ( 0 ) = 0 dt Mp ; dp k dt' 3 P j k £ x 3 z d p d P j 79 The gradient of F(p) i s then calculated as % . K n dUnx (t )) * ^ " - J i i o I ( 1 , k > 2 ( t a ^ ' t o i k ) -pj Integration of the state eguations was by i m p l i c i t integration using the second order trapezoidal rule}" 45] with future solutions obtained by a Newton-Raphson i t e r a t i o n . F u l l advantage was taken of sparsity i n the Jacobian of state equations with respect to s t a t e s i n simultaneous l i n e a r algebraic equation solutions within the Newton-Raphson i t e r a t i v e s o l u t i o n . Use of a l i n e a r inteqration method allowed simple coupling of the solution of the nonlinear state equations to the l i n e a r state s e n s i t i v i t y equations required to c a l c u l a t e gradients. In i n i t i a l diaqonal p o s i t i v e d e f i n i t e Hessian approximation qiven by . K n rd(£nx (t ))\2 hjj = h I I Ki.k) 2\ ^ k 1 3.7 d(.& .< - I i . i o 2 - l i t w k=l i=0 1 P J was supplied as required by Variable Metric methods, and t h i s provided much better s t a r t i n g c h a r a c t e r i s t i c s than the default i n i t i a l i d e n t i t y matrix estimate hjj =1. The Hessian approximation used i n Gauss-Newton l e a s t squares algorithms was not found e f f e c t i v e as a s t a r t i n g metric since i t was often singular or so badly conditioned as to be e f f e c t i v e l y so and therefore not s u i t a b l e . 3.1.5 I n i t i a l Parameter Guess And Obtaining A Good Starting Point The i n i t i a l parameter guess was to be obtained by assuming 80 a l l environmental functions constant at a reasonable mean l e v e l as estimated from a p r i o r i knowledge and gross inspection of nutrient and consumer l e v e l s . Since t r i a l parameter values s p e c i f i e d by the minimization procedure were unconstrained and could r e s u l t i n parameters y i e l d i n g unbounded solutions, p a r t i c u l a r l y i n the f i r s t few i t e r a t i o n s , some means of evaluating the solutions of these i n i t i a l wild gropings was reguired. This was achieved by a modified solution x(t,v,p,X), where X i s the set of observations x i k k=1,K i=0,n. The modified s o l u t i o n replaces the soluti o n x(t,v,p), the solution of a si n g l e i n i t i a l value problem, by the solu t i o n of a se r i e s of i n i t i a l value problems. This i s a s i m p l i f i c a t i o n of the s t a r t i n g procedure of van Domselaar and Hemker[47], At each observation time the model biomasses were set to observed values, thereby r e s t r a i n i n g otherwise unbounded solutions. This was a minor modification to the procedure of 3.1.4. This objective function F(p) was minimized to low accuracy using the same minimization procedure and the r e s u l t i n g parameter values used as an i n i t i a l point for s t a r t i n g the minimization of F{p). 3 T1.6 Termination Of Minimization Normal termination of the minimization procedure was on the condition that a l l parameters of each environmental function be changing between i t e r a t i o n s of the minimization procedure by l e s s than a s p e c i f i e d (EPSOPT) proportion of the i n i t i a l assumed mean value of the environmental function. After termination the 81 procedure was restarted i f the norm of the gradient exceeded a s p e c i f i e d value. This provided some check against premature termination. Limits on number of i t e r a t i o n s of the minimization were also provided. The program FMIN provided intermediate s t a t i s t i c s on the operation so that the progress of an minimization could be monitored and interrupted from an i n t e r a c t i v e terminal or recorded with batch output. Time l i m i t a t i o n s were not e x p l i c i t l y included i n the program controls, but were imposed with the operating system of the computer used. 3.1.7 Evaluation Of Parameter Estimates: Graphical Display Although the value of the r e s i d u a l function F (p) at termination of the optimizer provided a s i n g l e number describing the goodness of f i t of the r e s u l t i n g model, other c r i t e r i a were to be considered. Graphical display of the f i t t e d model solu t i o n , biomass data points, environmental smoothed input, and environmental data points was used to allow more judgement to be applied i n evaluating the r e s u l t i n g f i t . This was necessary to detect systematic departures of model s o l u t i o n from observations or e r r a t i c behavior of the model solution between data points. Graphical display of the polynomial description of estimated environmental functions was also used to evaluate the adequacy of estimation i n a more d i g e s t i b l e format than polynomial c o e f f i c i e n t s . D i f f e r e n t i a l biomass equations s o l u t i o n for the display was by fourth order Bunqe-Kutta inteqration (OBC BKC) error 82 co n t r o l l e d to 0.00001 r e l a t i v e error. Thus an independent check was a v a i l a b l e on the adequacy of the inteqration method used on the minimization procedures. 3.2 I d e n t i f i a b i l i t y Of Model: Noiseless Observations 3.2.1 Test Data Fiqure 3.1A displays the t e s t system, and Figure 3.1B the test environment, response of the te s t system, and observations to be used to te s t the i d e n t i f i c a t i o n procedure. The test system i s the f i r s t system of Example 3, Chapter I I . I n i t i a l conditions were taken from the steady state periodic solution of t h i s example. The s i n u s o i d a l temperature v a r i a t i o n has been perturbed to excite the nutrient so that i t i s not always close to a quasi-equilibrium of nutrient K(x,v). Samplinq i n t e r v a l s are 14 days for consumers, and 7 days for nutrient and environment with bursts of more rapid samplinq at 1 day i n t e r v a l s . There are K=65 biomass observation times. 3.2.2 I n i t i a l Guess And Assumed Structure The i n i t i a l parameters supplied to the minimization procedure i n t h i s and following cases consisted of a correct zeroth c o e f f i c i e n t i n a l l polynomials and zero c o e f f i c i e n t s for a l l hiqher terms. The zeroth c o e f f i c i e n t of the polynomial expansion i s the mean value of the function i n the Tchbycheff scaled ranqe, so t h i s i s equivalent to assuminq an i n i t i a l function constant at the mean value of the co r r e c t function. The number of unknown hiqher terms i n the polynomial expansion of each function was qiven as i d e n t i c a l to the number of terms 83 84 0.4 10 -1 10 -2 10 -4 NUTRIENT 0.0 103.0 ilflTo 1 U . 0 220.0 2CO.0 300.0 340.0 MO.O LOG SPECIES COLD IOO.O i<o.o teo.o 220.0 260.0 ORIS 300.0 140.0 3 M . 0 io - 1H 10 -2 30 ENVIRONMENT LOG SPECIES WARM 100.0 141.0 163.0 220.0 260.0 DATS 300.0 340.0 J M . O FIGURE 3.IB: Response of Test System to a particular environmental signal. ENVIRONMENT shows environmental observations and the spline function used as the environmental signal. NUTRIENT, LOG SPECIES COLD, and LOG SPECIES WARM time responses are shown along with sampled observations to be used i n estimating parameters. 85 i n the true function. This gave 16 unknown parameters of which 8 had correct i n i t i a l values. 3.2.3 Minimization Performance Table II gives the performance of three stages of minimization. The i n i t i a l approximation using F (p), a closer approximation c a l c u l a t i n g F (p) with i n t e g r a t i o n step egual to the i n t e r v a l between successive biomass observation times, and the f i n a l determination c a l c u l a t i n g F(p) with an inteqration step of .25 days. The convergence rate of the f i n a l minimization was greater than l i n e a r , the error function F(p) being reduced by close to one order of magnitude every n=16 i t e r a t i o n s . The estimated computation time f o r one integration time step was 0.00165 sec to c a l c u l a t e 3 new states, 48 new state s e n s i t i v i t i e s , and 16 gradient increments. The estimated computation time per i t e r a t i o n of the minimization exclusive of F(p) evaluation was 0. 1628 sec. Parameter values were determined within four f i g u r e accuracy i n the f i t . Other runs from d i f f e r i n g s t a r t i n q positions and using a d i f f e r e n t objective function also produced correct parameter values, e s t a b l i s h i n q the i d e n t i f l a b i l i t y of the model with these test data. 3.3 Estimation Of Model From Noisy Biomass Observations The biomass observations of the previous example were multi p l i e d by a loqnormally d i s t r i b u t e d pseudo-random variable 5=exp(n) where n was normally d i s t r i b u t e d with mean zero and standard deviation .35. The random variables were generated by 86 TABLE I I : Optimizer S t a t i s t i c s : Exact Observations F i r s t F i n a l Max Est I t e r s Integ Comp EPSOPT Term Value Value Integ No Time F(P) F(P) Step Integ Sec rmsre rmsre Size Days Steps 0.3 & 1.73 0.0065 & 1.08 7 65 9 11 2.77 10-1 A * 0.032 1.20 0.0036 1.06 7 65 46 79 16. 10-i* L * 0.022 1.15 .15x10-* 1.00039 .25 984 107 108 193. 10-s A F i r s t Value - Deviations of solution for i n i t i a l parameter estimate as F(p), mean sguare deviation of biomass natural logarithms and root mean square r a t i o error. F i n a l Value - Deviations of f i n a l estimate Max Integ Step Size - Time step of i n t e g r a t i o n used provided observation not encountered. Est No Integ Steps - Number of integration time steps reguired to calculate value of F(p). This i s a function of integration step s i z e , t o t a l time i n t e r v a l , number of observation times, and rate of change of the so l u t i o n . I t e r s - Number of i t e r a t i o n s of FHIN u n t i l termination Integ- Number of evaluations of F (p) reguiring one integration u n t i l termination. Comp Time- Execution time of procedure coded i n FORTRAN H . on IBM dual 370/168 machine. EPSOPT - Reguested accuracy of estimate, see 3.1.3 for d e t a i l s . Term - Condition upon which minimization was terminated: A - requested accuracy attained. L - maximum integrations exceeded T - allocated time exceeded S F^p) i n i t i a l function * Differences i n f i n a l and s t a r t i n g values for the same parameter estimate r e f l e c j ^ differences i n d i f f e r e n t i a l equation solution between F'(p) and F (p) with d i f f e r e n t i n t e q r a t i o n stepsize. NUTRIENT ENVIRONMENT OfllS -1 LOG. SPECIES COLD 10 -1 10 -2 IO9.0 l«.0 IM.O 120.0 KD.t 300.0 J«.0 MO.O OfltS i 10 -4 LOG SPECIES WARM 100.0 149.0 IM.O JM.0 2C0.0 OATS Jao.o 3«.o FIGURE 3-2: Noisy rapid observations and Test System response. The continuous biomass solutions NUTRIENT, COLD, and WARM are the response of Test System to the input signal ENVIRONMENT. Observations shown on the biomass graphs are the true biomass values multiplied by a lognormally distributed psuedo-random variable $ = exp(^) where ^ i s normally distributed with mean zero and standard deviation 0.35 • Nutrient was sampled at some times with a frequency of one per day. 88 UBC RANDN. Figure 3.2 displays the true s o l u t i o n and noise corrupted observations. The expected value of F(p) f o r the true s o l u t i o n i s <n 2> = .1225 or 1.419 rmsre. Table III gives the performance of a three stage minimization. TABLE I I I : Optimizer S t a t i s t i c s : Noisy Observations F i r s t F i n a l Max Est Iters Integ Corap EPSOPT Term Value Value Integ No Time F(p) F(P) Step Integ Sec rmsre rmsre Size Days Steps * 0.44 & 1.94 0.159 S 1.49 7 65 21 30 10-» A * 0. 15 1.49 0. 100 1.37 7 65 39 40 8.3 10-s L * 0. 126 1.43 0.92 1.35 .25 984 100 103 183. 10-s T £ F(p) i n i t i a l function * Difference i n f i n a l and star t i n g values for the same parameter estimate r e f l e c t ^differences i n d i f f e r e n t i a l eguation solution between F(p) and F(p) with d i f f e r e n t i n t e g r a t i o n stepsize. Figure 3.3B displays the r e s u l t i n g s o l u t i o n of the f i n a l parameter estimate., As evident, the value of F (p) attained by the f i t t e d solution i s lower than the expected value of F (p) for the true s o l u t i o n , i n d i c a t i n g that parameters may be underdetermined. Comparison with the true s o l u t i o n i n Figure 3.2 reveals that the l e v e l s of species HARM between days 160 and 200 have been underestimated and the early l e v e l of COLD has been underestimated between 160 and 180. This i s guite reasonable i n l i g h t of the perturbed observations. From Figure 3.2 i t i s quite evident that a run of noisy observations of WARM i n t h i s time i n t e r v a l have been by chance a l l below the true 89 FIGURE 3-3A: System function estimates from noisy rapid observations. 10.9 0.0 NUTRIENT o U3.Q 140.0 163-0 333.0 0_qfg° 343.0 323.3 90 211 ENVIRONMENT 100.0 14).0 3 0 0 . 0 3 C I . 0 3 4 3 . 0 10 -1 10 -2 10 -4 LOG SPECIES COLD o IC3.0 140.0 1S3.0 arcs :x>.s X3.3 340.0 360.0 10 -1. 10 -2 10 -4 LOG SPECIES WARM o 1 1 1 1 1 r 103.0 140.0 . 103.0 730.0 « 0 333.0 3C CflfS FIGURE 3.3B: Biomass response of Figure 3«3A estimated system to ENVIRONMENT. Observations of biomass shown are the noisy rapid observations used i n the f i t . 91 model s o l u t i o n . Errors i n nutrient estimate are p a r t i c u l a r l y noticeable i n the sharp spring peaks from day 120 to day 200, but errors are noticeable throughout the solution period. These errors are understandable since the rapid nutrient changes evident i n Figure 3.2 are not reasonably represented where only one week samples are taken. Notwithstanding these errors, the attained estimate of the solution i s remarkably representative of the true so l u t i o n and i s believed to be f a r superior to the r e s u l t attainable by any smoothing method which does not make use of the (known) d i f f e r e n t i a l eguation structure, environmental dependence, and implied continuity of the sampled variables. Now consider the accuracy of the estimated functions i n Figure 3.3A. Although e r r o r s i n magnitude and shape are evident i n the estimated functions t h e i r resemblance to the o r i g i n a l functions (Figure 3. 1A) i s c l e a r . The maximum error of composite functions r ( v ) / b ( v ) , e^ (v) a i (v) a^ (v)/b (v) are smaller than the errors of t h e i r component functions taken separately. The function m(v)tfARH has been evaluated to greater than graphically distinguishable accuracy, and m(v)COLD has been evaluated within a maximum error of approximately one t h i r d of the correct average value. The g i s t of t h i s b r i e f analysis i s that the composite functions comprising the memoryless nutrient approximation to the eguations have been more adequately estimated than the parameters describing rapid nutrient behavior. 92 Errors i n a(v) functions are most pronounced at low values of v. This may be attributed to the lesser number of samples at low values of v, but more c r i t i c a l seems the f a c t that low values of v occur when consumers have low values. Their e f f e c t on the nutrient would also be small at these low values of v, and therefore t h e i r influence i s masked by observation noise at low v, r e s u l t i n g i n poor d e f i n i t i o n of a(v) i n that range. Generally large errors i n a(v) are also to be expected since i t i s a c o e f f i c i e n t in the nutrient equation which by previous considerations may have underdetermined components related to the f a s t nutrient dynamics. The general conclusions of t h i s t e s t with error prone observations are that errors in the functions involved in the nutrient eguation are to be expected due to the p o t e n t i a l rapid dynamics of that variable r e l a t i v e to sampling rate, and that small e f f e c t s of consumers on nutrient at low temperatures are obscured by observation errors. Since such detailed knowledge as was made use of i n t h i s example of both the "true" continuous biomasses and environmental functions w i l l be unavailable f o r natural systems, a prediction for a s i m i l a r s i n u s o i d a l environmental s i g n a l d i f f e r i n g i n detailed f l u c t u a t i o n s w i l l be made. A measureable c r i t e r i o n f o r the adeguacy of parameter estimates which does not r e s t on unobservables i s therefore definable i n terms of the p r e d i c t i v e c a p a b i l i t y of the model. Figure 3.4 displays the new environmental s i g n a l and solution of the f i t t e d model equations. 93 0.72-, 0.0 LOG SPECIES COLD LOG SPECIES WARM FIGURE 3 . 4 : Predicted response of the estimated system of 3-3A to differing environment. Observation points are sampled from the true response of the Test System to the displayed environment. Data points shown i n Figure 3.4 are exact observations of the response of the true test system to the new environmental s i g n a l . Some errors i n the nutrient sol u t i o n are present, but consumer biomasses are well described despite errors i n nutrient and nutrient c o e f f i c i e n t s . I t i s granted that there w i l l be severe prediction errors i f the environmental signal i s very d i f f e r e n t from the one used i n f i t t i n g the model. However in f i e l d s i t u a t i o n s i t i s to be expected that there w i l l be strong deterministic components of the environmental s i g n a l and thus the model may be quite adeguate for prediction provided a representative sample of environmental inputs and biomass l e v e l s are available even though some parameter values have high e r r o r s . 3.4 Estimation Of Model From Slowly Sampled Observations: Indeterminate Parameters The rapid observations of section 3.2.1 taken on the nutrient at l e s s than 7 day i n t e r v a l s were deleted from the exact observations and an i n i t i a l approximation using F (p) made. Table IV gives the r e s u l t of the i n i t i a l minimization. The TABLE IV: Optimizer S t a t i s t i c s : Slow Observations F i r s t F i n a l Max Est I t e r s Integ Comp EPSOPT Term Value #alue Integ No Time F(P) F{p) Step Integ Sec rmsre rmsre Size Days Steps 0.42 0.083 7 34 15 18 3.13 10-2 A 1.91 1.33 next stage of the minimization was unable to improve upon t h i s 95 estimate and a high accuracy integration was not performed. Other s t a r t i n g points led to d i f f e r i n g parameter estimates with very s i m i l a r biomass solutions but much d i f f e r i n g function estimates primarily d i f f e r i n g i n the r(v) with g (v) , w(v), and u (v) composite functions well determined. Figure 3.5B displays the r e s u l t i n g f i t . Errors i n the functions (Figure 3.5A) are most apparent at values of v below 11. In order to eliminate indeterminacy of parameters and to v e r i f y the source of indeterminacy, the parameters of r(v) were fi x e d at the correct values and the minimization performed with the remaining free parameters. Table V gives the optimizer's TABLE V: Optimizer S t a t i s t i c s : Slow Observations, Fixed r (v) F i r s t F i n a l Max Est I t e r s Integ Comp EPSOPT Term Value Value Integ No Time F(P) F(P) Step Integ Sec rmsre rmsre Size Days Steps 0.237 0.0042 7 34 21 23 3.7 10-2 A 1.62 1.07 performance with r(v) parameters f i x e d at t h e i r correct values. The prediction of states as shown i n Figure 3.6B i s s l i g h t l y better than that of Figure 3.5B, but the unknown functions (Figure 3.6A) have been determined to a high degree of accuracy which could no doubt be made exact with higher integration accuracy. 3.5 Estimation Of flodel From Slow And Noisy Observations: Given r (v) Figure 3.7B displays slow sampled noise perturbed 96 97 FIGURE 3-5B: Biomass response of 3«5A estimated system to ENVIRONMENT. Observations fitted to are no more rapid than one week intervals. 98 1.2i FIGURE 3.6A: System function estimates from slow observations given correct r(v). 9 NUTRIENT FIGURE 3.6B: Biomass response of 3.6A estimated system to ENVIRONMENT. Observations shown and fitted to are no more frequent than one per week. 100 l . C h 1.2i O . l n 6.Ch 0 . 1 , FIGURE 3.7A: System function estimates from slow noisy observations given correct r(v). 101 1.0 0.0 NUTRIENT 4 03 o o 1 1 i 1 1 1 i tt-3.0 143.0 1 W . 0 273 .0 3*0.0 333.0 343.0 3a3.0 DRTS 21, 134 ENVIRONMENT 1 1 1 ! "I 1 100.0 143.0 183.0 320.0 2S3.0 333.0 343.3 3 . ^ DnlS LOG SPECIES COLD 10 -1 10 -2 10 -4 - i r iiia.o t«o.o i».o 730.0 rf.: a 330.0 KJ.O 300.0 ORIS 10 -1 10 -a 10 -4 LOG SPECIES WARM — i 1 1 1 1 1 — 100.0 140.0 103.0 270.0 7S3.0 333.0 340.0 0SI5 FIGURE 3.7B: Biomass response of 3-7A estimated system to ENVIRONMENT. Observation points shown and fitted to are response values of the true system multiplied by a lognormally distributed psuedo-random variable <f = exp(f) where ^ is normally distributed with zero mean and standard deviation .35 • Observations are no more frequent than one per week. 102 observations of the test system and the r e s u l t i n q biomass f i t . Since r(v) was badly determined without noise, f i t t i n g was performed using the r(v) parameters f i x e d at the correct values. Table VI displays the optimizer performance. Bote that the f i n a l F (p) i s l e s s than the expected value .1225 TABLE VI: Optimizer S t a t i s t i c s : Slow Noisy Observations, Fixed r (v) F i r s t F i n a l Max Est Iters Integ Comp EPSOPT Term Value Value Integ No Time F(P) F(P) Step Integ Sec rmsre rmsre Size Days Steps 0.49 0.20 7 34 9 14 2.7 10-2 & 2.0 1.56 0.405 0. 1201 7 72 42 60 11.8 10-s L 1.89 1.41 0.1372 0.1109 1 246 49 50 24.4 10-s x. 1.45 1.39 fo r the correct parameters. 3.6 Results The unconstrained minimization procedure chosen was e f f e c t i v e i n obtaining a minimum of the objective function of 16 parameters in a l l cases. The objective function F(p) was not observed i n any test cases to have d i s t i n c t l o c a l optima leading to very poor f i t s . Convergence to d i f f e r e n t parameter estimates with nearly equivalent biomass solutions was assumed to be a sign of underdetermined parameters. Parameters were underdetermined when observation noise or i n s u f f i c i e n t l y rapid nutrient sampling allowed large e r r o r s to be made i n c o e f f i c i e n t s r(v) and a i(v) with conseguent i n t e r a c t i o n with other parameters. Correlation of low v l e v e l s with low biomass 103 l e v e l s make functions inaccurate at low v l e v e l s . 3.7 Discussion The number of terms i n each function which may be s i g n i f i c a n t l y estimated i s a s t a t i s t i c a l guestion outside the scope of t h i s work, but the determinacy of the parameters considered only as minimum points of the function F (p) was addressed by considering the uniqueness of minima, and appropriateness of the state s o l u t i o n . The test data were, for the slow sampled noisy test case, i n s u f f i c i e n t to specify a F (p) with a unique minimum due to indeterminacy of functions r<v), a(v). Since t h i s s i t u a t i o n i s probably s i m i l a r to the actual data, a more riqourous approach to confidence estimates on the c o e f f i c i e n t functions was deferred u n t i l either a model better matched to t h i s data set i s developed or data better matched to t h i s model becomes a v a i l a b l e . Under these circumstances, the approach of Hemker and van Domselaarf47] may be applied usinq the Gauss-Newton Hessian approximation, Quasi-Newton Hessian approximation or a n a l y t i c a l l y calculated Hessian to define a paraboloid approximation to F(p) at the optimum and so obtain confidence estimates. Hore elaborately, a post-optimal analysis about the optimium locating points on given contours of F(p) with certain assumptions may be used to define sets of parameters of egual l i k e l i h o o d since F{p) may be interpreted as a l i k e l i h o o d function. This set can then be used to construct graphical outer confidence l i m i t s on the c o e f f i c e n t functions provided values of F(p) above the optimium can be linked with a 104 p r o b a b i l i t y of deviation. Computational and programming e f f o r t f o r such an analysis i s comparable to that reguired for the minimization procedure. A small number of t r i a l s on the noisy t e s t data where terms were given much i n excess of the number of non-zero parameters produced the following behavior. F i r s t , a more rapid convergence to a minimum was obtained than in cases with less parameters. Second, non-unigue parameter estimates were obtained. Third, and most obvious, the biomass solutions became very wiggly, fl u c t u a t i n g strongly between observations but giving a low value of F(p). A zero r e s i d u a l was not achieved i n these cases since the d i f f e r e n t i a l eguations remain a constraint on the attainable state solution even when parameters are badly determined. 3.8 Conclusions The data for the planned application of the estimation method, i f close to the rates and noise assumed in the test system do not contain s u f f i c i e n t information about r(v) and rapid nutrient dynamics. Fixing r(v) at some function may be e f f e c t i v e i n producing unique estimates of parameters. The i n s u f f i c i e n c y of data to determine uniquely the parameters of the qiven model makes development of confidence estimates premature inasmuch as the analysis w i l l be l a r g e l y inapplicable to the f i e l d data. 105 CHAPTER IV FITTING HODEL TO FIELD DATA The previous chapter addressed q u a l i t a t i v e l y the adequacy of parameter estimation when the observed system actually had the assumed model structure. When the observed system may not have the assumed structure i t i s misleading to speak of parameter estimation unless the assumed structure adequately describes the system behavior. Before considering how qood a parameter estimate has been obtained, i t must be established that parameter values e x i s t such that the assumed structure can model the system behavior. The optimization procedure, to the extent that i t gives the best f i t t i n g choice of parameters as defined by the minimum of the objective function F(p), selects the values of parameters under which the assumed structure best models the observed system. I t may be viewed as a Procrustean approach which sel e c t s parameter values without constraint to make the assumed structure f i t observations according to the s p e c i f i e d objective function. However, on the basis of the experience of Chapter I I I that the optimization procedure a c t u a l l y f i n d s the set of parameters r e s u l t i n g i n a global minimum of F (p), deviations of the model from observations i n d i c a t e those features of the actual system which are outside the assumed model structure. I f those deviations are systematic, then the assumed model structure must be rejected or replaced by an alte r n a t i v e structure. Refinements of parameter estimation accuracy are 106 questionable u n t i l the model structure i t s e l f has been validated by production of a f i t where model deviations from data are non-systematic and consistent with a random perturbation to observations. This c r i t e r i o n of non-autocorrelated random errors i n state estimates as a measure of goodness of f i t derives from the assumption of uncorrelated m u l t i p l i c a t i v e lognormal observation noise on biomasses. Under the assumption of egual variance of noise, the minimization of F{p) corresponds to the maximum l i k e l i h o o d estimate of parameters and states. If the number of measurements i s large compared to the number of unknown parameters, the minimum of F(p) i s unigue, and the r e s i d u a l deviations of the logarithm of biomass observations from the logarithm of biomass estimated are non-autocorrelated, then under the assumptions, the parameter and state estimates are unbiased and consistent. Thus the presence of correlated deviations of biomass f a l s i f i e s the s t r u c t u r a l assumptions. Addition of more parameters corresponding to s t r u c t u r a l change may remove correlated errors, but when the number of observations i s not much much greater than the number of parameters, as w i l l be the case here, each a d d i t i o n a l parameter increases the bias i n estimates, and when parameter numbers become equal to number of observations zero residuals may be attainable. H ithin the general model structure of Equations 1.1 - 1.9, the necessity of choosinq the number of parameters used to approximate each environmental function creates a vast number 107 of substructures. The optimization procedure can, for a given d i s p o s i t i o n of parameters obtain the best values of parameters but the question of how many parameters to use and their a l l o c a t i o n among functions cannot be decided only on the basis of the minimized value of F (p). A p r i o r i choices of number of parameters i n each environmental function were made considering the l i m i t a t i o n s on determining parameters from the number of observations on consumer and nutrient biomass, and the requirement that attack functions a(v) be d i s t i n c t i v e so that the i n t e r a c t i o n matrix (Equation 1.43 f f ) be nonsinqular and coexistence possible. Models f i t t e d from d i f f e r i n q sub-structures were judged inadmissable i f c o e f f i c i e n t functions were non-positive since these models may produce unbounded biomass solutions. Choice among the admissable models was made by the subjective assessment of the degree of systematic departure from observations, the smoothness of solutions between data points, and f i n a l l y by the value of F(p). Although the a p p l i c a t i o n of the method to the considered model did not include estimation of i n i t i a l conditions of states since i n i t i a l observations were av a i l a b l e , estimation of i n i t i a l conditions may also he performed by i n c l u s i o n of them as a d d i t i o n a l parameters for which s e n s i t i v i t y d i f f e r e n t i a l equations can be calculated as f o r other parameters. The data were f i t t e d with eleven assumed substructures and one was selected by.the previous c r i t e r i a . The deviations of t h i s model are discussed and the parameter estimates are 108 > • compared to some independently obtained components of the environmental functions. Another selected model i s used to predict, given temperature data from a d i f f e r e n t year, the phytoplankton and consumer biomasses. Structural changes i n the general model for estimation are suggested on the basis of the f i t t i n g experience. 4.1 Data The consumer community considered was the f i l t e r feeding zooplankton of Eunice Lake, a small oligotrophic temperate lake described by Federenkof 48 ]. From sampling described by Buckingham[33], estimates of population densities and mean lengths of zooplankters at standard depths to 30 m sampled at two week i n t e r v a l s from Hay 2, 1974 to December 12, 1974 were ava i l a b l e . Hater temperature, ash free seston (suspended pa r t i c u l a t e matter) density and phytoplankton volume density were also available over the same period at standard depths to 6 m at a one week sampling i n t e r v a l . Biomass (dry weight) and length measurements were avail a b l e (H.C.Swift unpublished data) from samples of zooplankton taken i n 1975 from Eunice lake and neighboring P l a c i d , Gwendoline, and Katherine lakes i n the U.B.C. Research Forest. These measurements were used to define a weight-length r e l a t i o n f or the four species Daphnia rosea, Holopedium gibberurn, Diaptomus kenai, and Diaptomus t y r e H i . Figure 4.1 shows the weight-length measurements for Daphnia rosea. Where an i n s u f f i c i e n t range of lengths was a v a i l a b l e , a r e l a t i o n weight=k(length) 3 was DRY WEIGHT mm Figure 4.1: Daphnia neighbo: . rosea weight-length measurements from four 'ring lakes and assumed weight-length relation. 110 assumed. Figure 4.2 shows how t h i s was applied to obtain the weight-length r e l a t i o n for Holopedium gibberum. Biomass density i n units of mg dry weight/m2 of lake were calculated from the depth d i s t r i b u t i o n of consumer populations with lengths and the weight-length r e l a t i o n for each species. The variations i n observed weight for given length were v i s u a l l y estimated to correspond to a standard deviation producing an upper l i m i t twice the lower l i m i t of weight. This corresponds to a rmsre of 1.42 or a natural logarithm standard deviation of 0.35. The remainder of zooplankton biomass was comparable to the l e v e l of D._ t y r e l l i . Phytoplankton biomass density i n mg dry weight /m2 was calculated by summing volume density estimates over 6 m depth and assuming a conversion of 1/6 dry weight to weight of water of the given volume. The environmental variable was chosen to be the average of lake temperature over the top s i x meters. This choice provided a r e l a t i v e l y smooth environmental s i q n a l at the samplinq i n t e r v a l of one week. Six meters was the approximate location of the thermocline during s t r a t i f i c a t i o n and i t was f e l t that t h i s region would encompass most of the photosynthetic production due to the drop i n temperature at lower depths and rapid l i g h t e x t i n c t i o n i n the brown-stained water of Eunice lake. Temperature data for 1975 were avail a b l e f o r Eunice lake so I l l Figure 4.2: Holopedium gibberum weight-length measurements from three^ neighboring lakes and assumed weight-length relation w=kl. 112 that 1975 biomass l e v e l s could be predicted with a qiven model, but s i z e and population data were not a v a i l a b l e u n t i l a f t e r the model was f i t t e d to 1974 data. Thus the p r e d i c t i o n of 1975 biomasses was t o t a l l y independent of 1975 biomass observations and i s a f a i r t e s t of the p r e d i c t i v e power of a model despite the subjective c r i t e r i a used in s e l e c t i n q the "best" substructure. 4.2 F i t t e d Sub-structures Table VII summarizes the f i t t i n q r e s u l t s of 11 substructure choices and three further f i t s on a selected substructure. TABLE VII: Summary of Substructure F i t s F i t T o t l No. Pararas/Fn rmsre Eps T Fe Fns End Int FIG No. No. r b a e m opt Sec >0 Stp Param Size 1 42 2 1(F) 3 4 3 1.39 .001 84 104 N L 1.75 2 46 2 1(F) 3 4 4 1.69 .01 26 68 N A 7 — 3 31 2 1(F) 4 1 2 1.58 .01 32 96 Y A 7 4 .3 4 32 3 1(F) 4 1 2 1.59 .01 38 136 N A 7 — 5 34 5 1(F) 4 1 2 1.69 .01 15 49 N A 7 6 24 3(F) 1 3 1 1 1.56 .001 91 101 Y L 1 — 7 25 3(F) 2 3 1 1 2.01 .005 41 75 N L 2 — 8 32 3(F) 1 5 1 1 1.52 .01 62 59 Y A 1 — 9 40 3(F) 1 5 3 1 1.39 .001 126 103 Y L 1 4 .4 10 48 3(F) 1 5 5 1 1.35 .001 130 96 Y L 1 — 11 43 3(F) 4 5 3 1 1.37 .001 67 59 Y L 1 — 12 40 3(F) 1 5 3 1 1.43 .001 51 106 Y L 3.5 4 .5 13 40 3 1 5 3 1 1.34 10-* 136 180 Y L 1.75 4 .6 14 40 3 1 5 3 1 1.35 10-* 75 100 Y L 1.75 4 .7 Column 2 qives the t o t a l number of parameters i n the f i t . Columns 3 - 7 qive the number of parameters used to approximate the functions. The f i r s t f i v e structures allowed functions m(v) to vary 113 with lake temperature v. In order to obtain determined c o e f f i c i e n t s , the c o e f f i c i e n t of b(v) was f i x e d . These f i t s were distinguished by m(v) f o r most species decreasing at higher v and where any functions became negative, the m{v) functions fo r Holopedium and D^ . t y r e l l i were negative above 15 deg. Figure 4.3 displays the estimates and f i t t o data for the only admissable system of f i t s 1 - 5 . Structures 6 - 11 were f i t t e d with r (v) fixed at the same function. In substructure 6 only three parameters were allowed f o r a(v) functions. Similar shaped a (v) functions resulted with the summer drops i n Holopedium and D. t y r e l l i not followed by the estimates of biomass. Increasing the number of a(v) and e (v) parameters allowed lower rmsre to be obtained with d i s t i n c t i v e p o s i t i v e consumer functions. Substructure 9, Figure 4.4 was selected as the "best" model considering the value of F (p) attained, p o s i t i v e c o e f f i c i e n t functions, non-systematic consumer errors and an attempt to keep the number of parameters down. Neighboring structures 7,9,10 and 11 were not s u b s t a n t i a l l y d i f f e r e n t i n biomass estimates or c o e f f i c i e n t functions so that t h i s estimate i s not strongly s e n s i t i v e to s t r u c t u r a l v a r i a t i o n of a,e,and b functions. F i t s 12 and 13 were used to determine the e f f e c t of d i f f e r i n g r(v) functions on the function estimates and biomass estimates. Figure 4.5 displays the f i t obtained when r (v) i s assumed to have one quarter the magnitude assumed i n f i t s 6 -11. Figure 4.6 displays the f i t obtained when a l l parameters l/day 0.2 r(v) phytoplankton 114 0.0 2.0 b(v) phytoplankton m g-day 0.0 2.0 a(v) Holopedium 21.0 5.0 20 m g-day a(v) D. kenai a(v) D. t y r e l l i m g-day 0.0 e/e o (v) 20 n e/e o 20 e/e\ o 20 g/g >(v) 0 >(v) 21". 0 0.41 1 day 0.0 0.4i 1 day 0.0 0.4 1 day O.O m (v) m (v) 9 13 17 21 'C 5 9 13 17 21 C 5 9 13 17 21 °C v: Average temperature top 6m Figure 4.3A: Estimated ( f i t 3) Eunice lake phytoplankton and consumer functions of environmental variable v, average temperature of top 6 m. Function b(v) is fixed at 1.0. Estimated functions a(v) are general cubic polynomials (four free parameters). 115 64 mg/m2 Phytoplankton Biomass 32 21 17 13 9 5 v : Eunice iy74 Holopedium Biomass mg/m 100 10 D.kenai Biomass mg/m , Daphnia Biomass mg/m T D . t y r e l l i Biomass 100 10 100 1 10 100 180 260 340 day 100 180 260 340 day Figure 4-3B: Estimated ( f i t 3) 1974 Eunice lake phytoplankton and consumer biomasses with observations used i n f i t and environmental variable v(t). 0.4 t r(v) phytoplankton 1.6 D(V) phytoplankton 116 1/day m g-day" 0.25 1 day 0.0 0.0 0.251 day 0.0 m (v) 0.251 m(v) 1 day m (v) 0.251 m(v) 1 day 0.0 9 13 17 21 °C 5 9 13 17 21 °C v: Average temperature top 6 m 5 ~~9 13 17 21 °C Figure 4.4A: Estimated (f i t 9) Eunice lake phytoplankton and consumer functions of environmental variable v used for predictions. Function r(v) was fixed at the given function. Estimated functions a(v) are general quartic polynomials (five free parameters). Estimated functions e(v) are general quadratics. 117 100 180 260 340 day 100 180 260 340 day Figure 4.4B: Estimated ( f i t 9) 1974 Eunice lake phytoplankton and consumer biomasses with observations used i n f i t and environmental variable v ( t ) . r(v) 0 QQ . phytoplankton 0.321 b(v) phytoplankton 118 0.25 1 day 0.0 m ill 0.25-x m(v) 1 day 0.0 0.25 1 1 day 0.0 m (v) a(v) D . t y r e l l i 100 g/g 0.0 0 0.251 m(v) 1 day 0.0 5 9 13 17 21 "C 5 9 13 17 21 *C* 5 9 13 17 21 C Figure 4.5A: Estimated ( f i t 12) Eunice lake phytoplankton and consumer functions of environmental variable v. Function r(v) was fixed at the given function which i s one quarter the magnitude of the function assumed i n Figure 4.4A. Figure 4 .5B: Estimated (fit 12) 1974 Eunice lake phytoplankton and consumer biomasses with observations used in f i t and environmental variable v(t). 120 Figure 4.6A: Estimated (fit 13) Eunice lake phytoplankton and consumer functions of environmental variable v. A l l functions were estimated. „ 121 are free i n the f i t , and the s t a r t i n g estimate f o r optimization i s f i t 12, Figure 4.5. The f i t s 1 - 1 3 using F (p) gave egual weight to a l l observations, producing large errors i n phytoplankton estimates r e l a t i v e to observed phytoplankton range of observations. I f i t i s assumed that the variance of observation errors on the logarithm of phytoplankton are one f i f t h of the variance of observation errors on the logarithm of consumers, then to obtain a maximum l i k e l i h o o d estimate, the sguared deviations of phytoplankton should be weighted by 5. Figure 4.7 ( f i t 14) displays the f i t when the objective function weighs phytoplankton sguared deviations of logarithm by a factor of 5. The general features of these f i t s as representative of the assumed model, may be summarized as: 1. A tendency for m(v) functions to be low or even negative at high temperatures for some consumers. ( f i t s 1 - 5, Figure 4.3) 2. Values of e(v) greatly i n excess of 1., the maximum e f f i c i e n c y possible i f the flow eguations are taken as energy compartment eguations. 3. A systematic error in a l l the nutrient estimates which overestimated nutrient i n days 160-180 and underestimaed nutrient i n days 280-300. During both these periods the average lake temperature was approximately 13 deg. Otherwise, nutrient estimates were not p a r t i c u l a r l y good, nor was there consistency i n t h e i r pattern. 122 Figure 4.6B: Estimated (fit 13) 1974 Eunice lake phytoplankton and consumer biomasses with observations used on f i t and environmental variable v(t). 123 0.0 1.0 m g-day 0.0 1.0 2 m g-day 0.0 1.0 m g-day 0.0 phytoplankton a(v) Holopedium a(v) Daphnia a(v) D.tyrelli 0.32" I 2 m g-day 0.0-»(v) 160 160 g/g. »(v) 160 g/g 0.4| m(v) 1 day 0.0 0.4i m(v) _1 day] 0.0 0.4 1 day m (v) 0.0 0.4 i m(v) 1 day 0.0 5 9 13 17 21 C 5 9 13 17 21 °C 5 9 13 17 21 °C Figure 4.7A: Estimated (f i t 14) Eunice lake phytoplankton and consumer functions of environmental variable v. A l l functions were estimated. This estimate was produced when squared deviations of logarithm phytoplankton biomass were weighted by five. 124 64 2 mg/m' 32 Phytoplankton Biomass mg/m Holopedium Biomass ioo A 10 l J • i mg/m£ 100 A 10 l A Daphnia Biomass 21 17 -I 13 9 -I v: Eunice 1974 mg/m n D.kenai Biomass 100 10 i A mg/m • 100-J 10 D . t y r e l l i Biomass 100 180 260 340 day 100 180 ' 260 ' 340 day Figure 4-7B: Estimated ( f i t 14) 1974 Eunice lake phytoplankton and consumer biomasses used i n f i t and environmental variable v(t)^ These estimates were produced when squared deviations of logarithm phytoplankton biomass were weighted by f i v e . 125 4. Given s u f f i c i e n t parameters in a{v) functions, consumer biomass was f i t t e d without systematic e r r o r . 5. Indeterminacy of r(v) magnitude with conseguent indeterminacy within a sc a l i n g factor i n a, e, and b was manifest but shapes of functions remained r e l a t i v e l y constant under d i f f e r i n g substructural assumptions, and values of t h i s s c a l i n g factor. The composite function r(v)/b(v) was, for example of consistent values. 4.3 Comparison Of Functions To Independent Measurements Table VIII compares f i t t e d functions of f i t 14 (Figure 4.6) i n the temperature range 13 - 17 deg. corresponding to times 215-260 with independently obtained measurements f o r some of the species. The l e v e l f o r r(v) was chosen so that a (v) functions would be consistent with estimated f i l t e r i n g rates. From the table, i t i s concluded that attack functions a(v) are not inconsistent with observed i n s i t u f i l t e r i n g rates of the consumers, nor are the m(v) functions inconsistent with consumer biomass turnover rates given the generation times and predation i n t e n s i t i e s on the consumers. The e f f i c i e n c y values e(v) are simply unbelievable as assi m i l a t i o n e f f i c i e n c i e s , since eff i c e n c y of assimilation cannot exceed 1 by thermodynamic considerations. 4.4 Interpretation Of F i t s The most serious deviation of functions i s the unreasonably high value of e f f i c i e n c i e s estimated. Although the method used t o measure phytoplankton density may have s e r i o u s l y 126 TABLE VIII: Comparison of Estimated Functions to Measured Rates Species Estimated Mean a(v) v 13-17 m2/g-day Estimated F i l t e r i n g Rate a(v)*6 m ml/yg-day Measured F i l t e r i n g Rate water temp 14-15 C ml/y g-day Holopedium 0.16 gibberum Diaptomus 0.07 kenai Daphnia 0.15 rosea Diaptomus 0.45 t y r e l l i .96 .42 .9 2.7 1. 1 [33]* 0.6 [33]* 1.4 [33]* 0.4 [33]* Holopedium gibberum Daphnia rosea Diaptomus kenai Diaptomus t y r e l l i Estimated Rate m(v) 1/day 32 ,25 24 25 Predation Rate 1/day Biomass Turnover Rate 1/day } . 1 - .05 [50] } -28peak- .07ave.*[49] } .07 -.1 peak- .03ave.* } T49] 006 [50] * measured i n Eunice Lake underestimated the density of larger c e l l s (C.J. Halters), t h i s error i f corrected could not increase phytoplankton observations by more than 10 times which would s t i l l reguire that e f f i c i e n c i e s be of order 10. The p o s s i b i l i t y of the high e f f i c i e n c y estimate being an a r t i f a c t of the f i t t i n g procedure i s rejected f i r s t l y since the approximately correct values of m (v) and a(v) with observed values of phytoplankton reguire e f f i c i e n c i e s of t h i s magnitude to maintain the consumer biomasses and second since there i s no p o s s i b i l i t y of dynamic e f f e c t s being neglected since the 127 d i f f e r e n t i a l equations are s a t i s f i e d by the chosen f i t t i n g method. I t i s therefore concluded that the assumption that phytoplankton form the sole nutrient of the consumers must be rejected. A nutrient supply approximately 100 times as dense must be present i f the attack rates a(v) and consumer biomass turnover rates m(v) are accepted and a s s i m i l a t i o n e f f i c i e n c i e s l e s s than one are to be estimated. The concentration of ash free dry weight of seston i s approximately 100 times the phytoplankton biomass observations and therefore constitutes a l i k e l y source of the missing energy. This a d d i t i o n a l food source might also account for the tendency of m(v) to become negative at higher temperatures since a negative value of n (v) i n the assumed structure formally implies an energy flow into consumers not related to phytoplankton density. I f t h i s food source were more important at higher temperatures, t h i s could account for the decrease of m(v) at higher temperatures. The unusual behavior of the f i t t e d a(v) functions at low temperatures i s probably due to the same factors which produced the same form of errors i n the test system f i t s i n Chapter I I I . The systematic errors of l a t e spring (days 160-180) and f a l l (days 280-300) phytoplankton estimates may be interpreted as an a r t i f a c t of the f i t t i n g procedure r e f l e c t i n g another fundamental v i o l a t i o n of the assumed structure i n these periods. Consider consumer growth i n these time periods. Consumer observations i n the l a t e spring a l l r e f l e c t very rapid growth 128 which has been f i t t e d by the model. Consumer observations i n the f a l l r e f l e c t declining biomasses which have again been f i t t e d by the models. Yet temperatures i n both periods are 12 -13 deg. So temperature induced changes cannot account through e f f e c t s on functions for the difference i n consumer growth i n these periods. Neither can declining phytoplankton food supply account for t h i s difference between spring and f a l l since the period of rapid spring growth has only three quarters of the phytoplankton density of the period of f a l l decline. Thus these observations are inconsistent with the assumption of a phytoplankton food supply whose e f f e c t s on growth are a function of temperature only. The f i t t i n g procedure, i n order to account for these differences i n consumer growth at the same temperature has produced increased nutrient density i n the sprinq period and decreased nutrient density i n the f a l l period so that consumer growth i s consistent with the assumed model structure. The obvious assumption of a seston food supply w i l l not correct these discrepancies i n the food density-temperature explanation of consumer growth i n those periods since the seston densities are also lower i n the sprinq than f a l l . The f a l l declines may be partly explained by predation. The estimates of Federenko[49 ] place the peak of predation by the major predators, the larva of Chaoborous amexicanus and Chaoborous t r i v i t t a t n s (Eunice lake had no f i s h p r i o r to 1975) at day 150 with a l e s s e r peak at 210-230 in years 1971 and 1972. B o r t a l i t y due to predation was maximum durinq the sprinq high 129 growth period and very low from day 240 u n t i l year end. The maximum predation estimated by Federenko was .28 per day (numerical standing crop) on Diaptomus kenai near day 150, 1972 and season average .11, with the other pery considered, Diaptomus t y r e l l i maximum predation of .1 per day (numerical standing crop) at day 150 with season averaqe of .03 per day. Predation was below average i n the f a l l period of growth decline for both prey and both predators. However, removal of the Chaoborous by introduced f i s h i n 1976 resulted i n increases i n the l e v e l s of Daphnia and kenai i n l a t e f a l l with a reduced rate of decline[53]. Noting that there may be changes i n the guality of the as yet poorly defined food source, t h i s p o s s i b i l i t y was not explored due to lack of evidence. The rapid spring growth cannot be explained by a so-called time delay i n response of the populations to the e a r l i e r high nutrient l e v e l s since the biomass measure of the consumer populations accounts immediately f o r any energy inflows i n contrast to a numerical measure of the consumer population which might show a delayed reproductive e f f e c t . The simplest v i o l a t i o n of the s t r u c t u r a l assumptions known to be present i s that a l l four of the consumers produce resting eggs which overwinter and hatch i n the spring, thus v i o l a t i n g the assumption of intraspecies biomass eguivalence and adding an unaccounted biomass flow and store. Although i t might be assumed that t h i s biomass contribution i n the spring may be 130 nutrient, the cost to the consumers of producing the resting egq biomass in the f a l l would be substantial and may account for the decline i n consumer biomass before a s s i m i l a t i o n l i m i t a t i o n might be i n f e r r e d from temperature and nutrient density. The s t r a t e g i c advantage of a n t i c i p a t i n g with some safety marqin, say one generation of 10-14 days the partly predictable f a l l food a v a i l a b i l i t y decline i s that biomass can be e f f e c t i v e l y t r a s f e r r e d to a safe store rather than expended i n maintaining a starving active form aft e r food becomes' unavailable. I t i s known[54] that production in Daphnia pulex of resting egqs i s increased by low temperature and shorter day length. A model continuous over a year period should therefore include a dormant l i f e stage for t h i s system. 4.5 Predictions Of The F i t t e d Model Neglecting the errors attributed to r e s t i n g egg production, the appearances of the basic nutrient-temperature s t r u c t u r a l hypothesis to explain consumer growth may be saved i f the nutrient i s assumed to be an unspecified food source which obeys the same dynamics as f i t t e d but whose units are 100 times the observed units. Figure 4.8 displays the predicted s o l u t i o n of f i t 9 (Figure 4.4A) to 1975 temperatures where i n i t i a l biomass values were taken from the same time in 1974. This substructure had been chosen before independent measurements and 1975 biomass data were avai l a b l e . , In gross terms, the summer temperatures i n 1975 were 131 200-1 Phytoplankton Biomass mg/m 100-1975 observations — model predicted biomass — 1974 observations v: Eunice 1975 0-mg/m 100.0 —! 1.0 0.01-Holopedium Biomass mg/m-i D.kenai Biomass mg/m^ 100.0 1.0 -0.01 • Daphnia Biomass , 2 + mg/m • * -«• . 100 D.tyrelli Biomass 200 day 300 i 100 200 day 300 Figure 4 .8 : Comparison of model biomass predictions (Figure 4 .4A, f i t 9) for 1975f given v(t), with 1975 biomass observations and 1974 biomass observations. The i n i t i a l biomasses were assumed at 1974 values. 132 >.• approximately 2 deq. cooler than i n 1974 with a s i m i l a r shaped but e a r l i e r temperature pattern. The phytoplankton model estimates are very much i n error past day 220, but since we have already rejected t h i s part of the structure, the error i s not s u r p r i s i n g . It i s i n t e r e s t i n g to note that the l e v e l of phytoplankton during the apparent bloom coincides with the quantity r ( v ) / b ( v ) , the l e v e l to which model phytoplankton would tend i f u n u t i l i z e d . The s i q n i f i c a n t differences predicted between 1975 and 1974 i n consumer biomass occur primarily a f t e r day 220 when the e f f e c t s of the lower summer temperatures on the model become apparent. Holopedium i s predicted to be hiqher by a factor of approximately 5 i n the period 220-300. Both Daphnia and D. kenai are predicted to beqin a decline around day 200 from t h e i r 1974 l e v e l s and to decline rapidly u n t i l they are reduced to .1 of 1974 l e v e l s . Predictions of the other f i t s i n Fiqures 4.3, 4.4, 4.5, and 4.6 were the same at t h i s l e v e l of precision of description of the consumers. Phytoplankton predictions were diverse but none showed the bloom apparent i n 1975 observations. He compare the effectiveness of the model as a predictor with simply the previous year's observations as shown i n Fiqure 4.8. The model prediction i s no better than the previous season i n predictinq Holopedium, Dy t y r e l l i , , and Daphnia v a r i a t i o n and i s worse i n prediction of D. kenai. The accuracy of prediction of f i n a l biomass i s worse fo r 133 D. kenai than the previous year. A prediction was also made where 1975 i n i t a l conditions were given for consumers. As Figure 4.9 shows, the prediction i s s l i g h t l y improved, but the same predic t i o n e rrors are present. Namely, an overestimation of the effectiveness of Holopedium at the lower summer temperatures with a connected drop i n Daphnia and D. kenai biomasses. This example serves i n c i d e n t a l l y to demonstrate the i n s e n s i t i v i t y of the model to i n i t i a l condition perturbations. Despite these errors i n detailed prediction, the model has p a r t i a l l y accounted f o r coexistence. The i n c o r r e c t l y predicted drops i n Daphnia and D. kenai are much les s severe than could be re a l i z e d i n a deterministic i n s t a b i l i t y . Table IX gives the solution of the f i t t e d model to three fixed environmental values with i n i t i a l conditions as given i n 1974 a f t e r 170 days. Three TABLE IX: Solutions of F i t 9 With Constant v. I n i t i a l Conditions - 1974 Time I n t e r v a l - 170 days F i n a l Biomasses v Holopedium D.kenai Daphnia D . t y r e l l i 17 .0001 50 400 * 0.1 16.25 8 .04 .08 160 * 14.4 160 * .004 .001 8 * only increasing species at end of time i n t e r v a l of the model species are demonstrated t o be capable of excluding others i n the respective fi x e d environments. The rates of exclusion of Daphnia and D. kenai f o r the two lower temperatures are much greater than the error i n growth f o r the predic t i o n of 200. mg/m2 10(H Phytoplankton Biomass 100 200 300 day © 1975 observations "~~ model predicted biomass 1974 observations v: Eunice 1975 100 .200 300 day m, g/m -i Holopedium Biomass 100.0-1.0-0.01 mg/m2-] 100.0 H 1.0 H o.or D.kenai Biomass mg/m' 10C.0 2 Daphnia Biomass l.Oi 0.01— 100 200 / 2 mg/m "I 100.01 1.0H 0.01-D.tyrelli Biomass 300 day 100 200 300 day Figure 4-9: Comparison of model biomass predictions (Figure 4.4A, f i t 9) for 1975 given v(t) with 1975 biomass observations and 1974 biomass observations. The i n i t i a l biomasses were assumed at 1975 values.. 135 . >.• maintaining the coexistence of the model at l e a s t and p a r t i a l l y accounting for the coexistence of the natural system. Periodic r e p e t i t i o n of the 1974 and 1975 temperature si g n a l s l e d , r e s p e c t i v e l y , to dominance of Daphnia and Holopedium. Finding periodic environmental si g n a l s to produce coexistent solutions of the f i t t e d models was not attempted i n view of the serious d e f i c i e n c i e s i n model accuracy and lack of continuity of measurements through winter. Repetition of 1974 temperatures did not lead to stable coexistence, but i t i s not cl e a r that the observed system would have a stable coexistence under t h i s regime since Daphnia over the period was by f a r the most successful in increasing biomass. 4 T6 Conclusions The model structure assumed fo r estimation was inappropriate to the data ava i l a b l e i n two ways. F i r s t , the u n s a t i s f i e d reguirement f o r rapid nutrient measurement introduced an indeterminacy i n function estimates due to slowly sampled phytoplankton and seston. Second, the necessity of having nutrient measurements when the i d e n t i t y of nutrient i s , i n t h i s case, problematical l i m i t s the a p p l i c a t i o n of estimation. The v i o l a t i o n of the equivalent species biomass assumption by r e s t i n q eqq production could i n p r i n c i p l e be avoided by r e s t r i c t i n q the application of the model to times before the f a l l decline i n temperature, but i f t h i s were done with the 136 a v a i l a b l e data, i n s u f f i c i e n t data would remain to define c o e f f i c i e n t functions i n any d e t a i l . These d i f f i c u l t i e s prevented a convincing t e s t of the extent or existence of environmental s t a b i l i z a t i o n i n the considered system. Nevertheless the use of the optimization procedure at the least demonstrated that a f i t can be achieved, and the discrepancies i n that f i t may be used to redefine the s t r u c t u r a l hypotheses without reference to p r e d i c t i v e f a i l u r e . The models predictions, while incorrect i n some respects were not b i o l o g i c a l l y absurd suggesting that modifications may yet be possible which would render the basic concept of an environmentally s t a b i l i z e d s i n g l e nutrient multiple consumer model testable and p o t e n t i a l l y applicable i n t h i s system. I t i s therefore suggested that the memoryless nutrient consumer system (Eguations 1.43) be used as a basis for further i n v e s t i g a t i o n . Since these eguations involve only consumer biomass and the environmental variable, nutrient observations are not reguired. The assumption of a guasi-eguilibrium of nutrient allows dynamics to be considered only at the consumer time scales which should speed up s o l u t i o n of the model d i f f e r e n t i a l eguations. Since there i s a degree of freedom i n the s p e c i f i c a t i o n of w and u function f o r the memoryless nutrient system to give i d e n t i c a l s o l u t i o n s , i t i s necessary to • remove t h i s before attempting i d e n t i f i c a t i o n . This i s done by choosing a reference species (species 0) and then defining an i d e n t i f i a b l e structure by making eith e r w (v) functions or u (v) 137 functions r e l a t i v e to the reference species. Eguation 4.1 and Eguations 4.2 give the memoryless nutrient eguations with the reference function w (v) before redefining r e l a t i v e functions. y = g (v) - -— T w u. exp(y.) Jo 6o w „ j ° 3 3 ° 3 w i y ± = g±(v) - — I w u exp(y ) 4.2 o j J J defining new functions which are r e l a t i v e to the defined u> = 1 o as w±(v) u i ( v ) = T O O 4 - 3 4.4 y ±(v) = W Q(V) u±(v) The i d e n t i f i a b l e form with r e l a t i v e functions i s 4 5 y = g(v) - ui(v) u'(v) exp(y_) where the reference species has u =1. o In performing these s i m p l i f i c a t i o n s , the only r e s t r i c t i o n s which remain from the o r i g i n a l r e s t r i c t i o n s on functions are that w (v) >0 and u(v) >0. The i d e n t i f i a b l e functions g, w and y have absorbed and confounded r e a d i l y i nterpretable b i o l o g i c a l processes to give an almost purely phenomenological description of consumer growth i n terms of the environmental variable and competing biomasses. Given the reference function and assumed values f o r some of the more basic functions, however, other basic functions may be recovered from t h i s s t r u c t u r e . This 138 formulation of model structure dictates a more riqourous s t a t i s t i c a l approach to i d e n t i f i c a t i o n since the v e r i f i c a t i o n of parameter values w i l l be dependent on these c r i t e r i a alone. With t h i s model, the data of Federenko from previous years 1971 and 1972 may be used although phytoplankton and seston were not sampled. Thus i t should be possible to t e s t the memoryless nutrient model over the times when r e s t i n g eggs may be neglected using three years 1971, 1972, and 1974 and obtain l e s s biased estimates of the functions. This i s a good d i r e c t i o n for further work. A summary of the conclusions of using the assumed model structure and f i t t i n g method are that consumers are not feeding primarily on phytoplankton. Temperature was shown to have s i g n i f i c a n t e f f e c t s on consumer growth, but errors i n s t r u c t u r a l assumptions and the small number of samples prevented a t e s t of the extent of i n t e r s p e c i f i c competition f o r a single nutrient with temperature s t a b i l i z e d coexistence. 139 CHAPTER V CONCLUSION Results Of Investigation A dynamic model has been given where multiple consumers of a single nutrient per s i s t i n a f l u c t u a t i n g environment. The model s a t i s f i e s e n e r g e t i c a l l y reasonable constraints on the component species, and does not allow unbounded biomass. Persistence i s produced by the e f f e c t s of a f l u c t u a t i n g environment on the a b i l i t y of each species to u t i l i z e and assimilate the single nutrient. Environmental fluctuations may be random or periodic to obtain persistent solutions. The persistence of examples resembling two natural systems was demonstrated. The form of the model allows i d e n t i f i c a t i o n of unknown species c h a r a c t e r i s t i c s which are functions of a measureable environmental variable given sampled noise free observations of biomass and the environmental variable. When non-autocorrelated m u l t i p l i c a t i v e lognormal noise of known d i s t r i b u t i o n corrupts biomass observations, the estimation method can produce reasonably accurate biomass predictions f o r an a r t i f i c i a l system. The model was tested using data from the grazer zooplankton community of a small lake. A f i t t e d model was produced which was bounded but was ene r g e t i c a l l y impossible and i t was concluded that phytoplankton cannot be the major food of the 140 ' >•-consumers. Despite t h i s , predictions of independent zooplankton biomass data were not absurd, but were l e s s accurate than assuming a r e p e t i t i o n of the previous y e a ^ s data. Comparison of model species persistence with observed species persistence was not possible due to the l i m i t e d span of data. I f large declines of model predicted biomass from observations are interpreted as exclusion i n the model, two species are excluded i n the prediction. However t h i s exclusion i s not so rapid as can be produced in the model, and coexistence of at l e a s t two species i s therefore accounted for by the f i t t e d model over the l i m i t e d time span. Suggestions For Further Work A proof of existence of a stable solution to the model should be sought. I t i s conjectured that every example of the model has a companion time invariant system which has i d e n t i c a l persistence properties for a l l 2 n species combinations and whose e q u i l i b r i a match the means of the model solutions. This companion system i s not equivalent to the r e l a t e d time i n v a r i a n t system used to obtain a mean approximation. Examples were qiven where the model had a persistent s o l u t i o n but not the related system, and the related system had a stable equilibrium but the model community was not persistent. I t i s also conjectured that the condition <C> nonsinqular i s a necessary condition for a persistent s o l u t i o n . This would c o n s t i t u t e a p r i n c i p l e of competitive exclusion for these systems. Examples for which the f e a s i b l e mean approximation i s 141 unstable should be constructed to determine whether decay of these systems i s i n i t i a l condition dependent. The i d e n t i f i c a t i o n technique should be used to attempt synthesis of models which approximate some qiven periodic solutions. Using t h i s technique some of the more i n t e r e s t i n g systems of r e a l i s t i c rates may be explored. As an example, a coexistent consumer community i n an almost periodic environment consisting of a reqular sequence of seasonally dominant environmental s p e c i a l i s t s with rare generalized species should be attempted. The question of int e r e s t in such a system i s whether a departure from the normal seasonal pattern could produce an i r r u p t i o n of a rare q e n e r a l i s t . The use of orthoqonal polynomials to estimate time varying co e f f i c e n t s of eguations should be investigated i n l i n e a r d i f f e r e n t i a l eguations. This y i e l d s a set of time invariant l i n e a r eguations with inputs T n(v) multiplying unknown but fix e d parameters. From t h i s , some guidance might be obtained for the number of terms estimable from sampled observations of a system. The i d e n t i f i a b l e form of eguations 4.5 should be applied to an available larger set of data from the U.B.C. Research Forest Lakes and a comparison made between the v a r i a b i l i t y accounted f o r by an average seasonal v a r i a t i o n , separate temperature dependent l o g i s t i c models, and temperature s t a b i l i z e d competition of the i d e n t i f i a b l e form. The s t a t i s t i c a l confidence i n a f i t t e d model f o r assumed noise s t a t i s t i c s may be 142 estimated by a Monte Carlo method. Taking the f i t t e d model as a te s t system, d i f f e r i n g noise patterns of equal p r o b a b i l i t y are generated and these patterns are used to corrupt the f i t t e d model observations to produce a set of randomly perturbed observation sequences. The v a r i a b i l i t y of parameter estimates from f i t s to each member of the set of observation sequences estimates the s e n s i t i v i t y of parameter estimates to random v a r i a t i o n , that i s , s t a t i s t i c a l confidence i n the model estimated. The v a r i a b i l i t y of biomass solutions among the various f i t s estimates the confidence i n biomass estimation, and may s i m i l a r i l y give a spread of predictions to estimate confidence in prediction given the assumed noise s t a t i s t i c s . The memoryless nutrient equations 1.43 also approximate the dynamics of a multiconsumer community i n a mixed, rapidly flowing chemostat. Thus the construction of persistent microcosms i n a con t r o l l e d periodic environment i s f e a s i b l e since there are quides to t h e i r construction and s t a b i l i z a t i o n . Such systems should be constructed to test the effectiveness of estimation procedures on a system f o r which many samples can be obtained. I f estimation proves e f f e c t i v e i n predictinq the e f f e c t s of planned manipulations, then such microcosms could be used to test r e v o l u t i o n a r y hypotheses i n a prec i s e l y describable and c o n t r o l l a b l e environment. The extension of the model to a s p a t i a l l y d i s t r i b u t e d system might be approached by assuming di s p e r s a l of consumers from an area as the product of density and a po s i t i v e siqmoidal 143 function of s p e c i f i c growth rate. However, .a simpler d i s t r i b u t e d system should f i r s t be devised before attempting the problem of dispersing consumers competing i n a s p a t i a l l y and temporally varying environment for a nutrient which i s dynamic. Competition for space i t s e l f i n a temporally f l u c t a t i n g environment should be formulated f i r s t with d i s p e r s a l of the assumed form. Competition between multiple l i f e stage orqanisms should be investigated f o r the s p e c i a l case when the e a r l i e s t l i f e stage of each species competes f o r a single nutrient i n a fl u c t u a t i n g environment. Such a model might explain the r e s u l t s of S u t c l i f f e et al.[52] who showed that variat i o n s i n sea surface temperature near the f i r s t growing season of each species explain at l e a s t 50% of the f l u c t u a t i o n in catch of many species, yet the t o t a l catch f o r 17 species over 19 years was within ±10% of the mean. These two r e s u l t s are consistent with the concept of. a change-stabilized consumer community lim i t e d by one resource. Conclusions Environmental f l u c t u a t i o n s should be considered as one of the bases f o r persistence of e c o l o g i c a l communities. The elemental model given i s one way of representing the persistence e f f e c t s of environmental f l u c t u a t i o n i n a t e s t a b l e manner. A conclusive example of a b i o l o g i c a l community s o l e l y change-stabilized as i s assumed i n the elemental model has not been given. The r e l a t i v e success of prediction with the model 144 i n application to grazer zooplankton of a small lake i s interpreted as evidence of p a r t i a l change-stabilization. 145 REFERENCES [1] Lotka, A.J., Growth of mixed populations, J. Hash. Acad. S c i . 22:461-469,(1932) [2] Lotka, A.J. Elements of physical biology, Baltimore: Williams and Wilkins. ( reissued as Elements of Mathematical Biology by Dover, 1956) (1925) [3] Volterra, V., Lecons sur l a theorie matheaatique de l a l u t t e pour l a 4 v i e , Paris: G a u t h i e r - V i l l a r s . (193 1) [4] Volterra, V., V a r i a z i o n i e f l u t t u a z i o n i del numero d'i n d i v i d u i in specie animali conviventi, Mem. Acad. L i n c e i . 2:31-133, (tr a n s l a t i o n i n an appendix to Chapmanfs Animal Ecology, New York, 1931) (1926) [ 5 ] Beverton, R.J.H., and Holt, S.J., On the dynamics of exploited f i s h populations, Min. Agr. Fish. And Food (U.K.), Fish. Investig. Ser. II (19):533, (1957) [6] Schaefer, M.B., A study of the dynamics of the fishery for yellowfin tuna i n the eastern t r o p i c a l P a c i f i c Ocean, Inter-Amer. Trop. Tuna Comm. B u l l . 2 (6): 245-285 (1957) [7] Verhulst, P.F., Notice sur l a l o i gue l a population s u i t das son accroissement, Corr. Math. Phys., Paris: A Quetetet. (1839) [8] Park, R.A., A generalized model f o r simulating lake ecosystems. Simulation 23(2):33 (1975) [9] Walters, C.J., and Peterman, R.W., A systems approach to the dynamics of spruce budworm i n New Brunswick, Quaestiones Entomologicae 10:177-186(1974) [10] Rosensweig, M. , Paradox of enrichment: der>tabilization of exploitation ecosystems i n e c o l o g i c a l time. Science 171:385-387 (1971) [11] May, R.M.,and MacArthur,R.H., Niche overlap as a function of environmental v a r i a b i l i t y , Proc. Nat. Acad. S c i . 69:1109-1113 (1972) [12] G i l p i n , M.E., Do hares eat lynx?. Am. (1973) Nat. 107:727-730 146 [13] Ayala, F.J., G i l p i n , M.E., and Ehrenfeld, J.G., Competition between species, t h e o r e t i c a l models and experimental t e s t s , Theor. Pop. B i o l . 4:331-356, (1973) [14] Hay, R.M., S t a b i l i t y and complexity i n model ecosystems, Princeton University Press, (1973) [15] Hay, R.M., S t a b i l i t y i n randomly f l u c t u a t i n g versus deterministic environments. Am. Nat. 107:621-650 (1973) [16] Feldman, H.W., and Roughqarden, J . , A population's stationary d i s t r i b u t i o n and chance of extinction in a stochastic environment with remarks on the theory of species packinq, Theor. Pop. B i o l . 7:197-207 (1975) [17] Long, G.E., Model s t a b i l i t y , r e s i l i e n c e , and management of an aquatic community, Oecologia (Ber 1. ) 17: 65-85 (1974) [18] Long, G.E., Duran, P.H., Jeffords, P.O., and Heldon, D.N., An application of the l o g i s t i c equation to the population dynamics of salt-marsh gastropods, Theor. Pop. B i o l . 5:450-459 (1974) [19] Maguire, B.,jr.,Ecosystem simulation through use of models of subsystem response structures. Simulation 23:(5):149-158, (1974) [20] Hutchinson, G.E., The paradox of the plankton. Am. Nat. 95:137-145, (1961) [21] Grenney, H.J., B e l l a , D. A., and Curl H.C., A t h e o r e t i c a l approach to i n t e r s p e c i f i c competition i n phytoplankton communities. Am. Nat. 107:405-425,(1973) [22] Stewart F.M., and Levin B.H., P a r t i t i o n i n g of resources and the outcome of i n t e r s p e c i f i c competition: a model and some general considerations. Am. Nat. 107:171-197, (1973) [23] Armstrong, R.A., and McGhee, R., Coexistence of species competing for shared resources, Theor. Pop. B i o l . 9:317-328, (1976) [24] McGhee R., and Armstrong, R.A., Some mathematical problems concerning the e c o l o g i c a l p r i n c i p l e of competitive exclusion, J . D i f f . Eq. 23:30-52 (1977) [25] Koch, A.L., Coexistence r e s u l t i n q from an al t e r n a t i o n of density dependent and density independent qrowth, J . Theor. B i o l . 44:373-386, (1974) 147 [26] Koch A.L., Competitive coexistence of two predators u t i l i z i n g the same prey under constant environmental conditions. J . Theor. B i o l . 44: 387-395 (1974) [27] S i n c l a i r , A.R.E, The resource l i m i t a t i o n of trophic l e v e l s i n t r o p i c a l grassland ecosystems, J . Anim. Ecol. 44:497-520, (1975) [28] J e f f r i e s , C., P r o b a b i l i s t i c l i m i t cycles, i n Mathematical problems in biology, V i c t o r i a conference 1973, [ed. Driessche, P. Van den]. Springer Verlanq 123-131, (1974) [29] Braun, H.M.H., Primary production i n the Serenqeti: purpose, methods, and some r e s u l t s of research, Annales de 1»Universite d'Abidjan, (E)6:171-188 (1973) [30] Wetzel, R.G., Limnology, W.B. Saunders, Philadelphia, 302,(1975) [31] i b i d . 578 [32] i b i d . 444,445 [33] Buckingham S. L., Ph.D. Thesis (unpublished) , University of B r i t i s h Columbia (1977) [34] Mayzaud, P., and Conover, R.J., Influence of pot e n t i a l food supply on the a c t i v i t y of digestive enzymes of n e r i t i c zooplankton, Proc. Of 10th European Symposium on Marine Biology, Ostend, Belgium, Sept. 1975, T (eds) Persoone, G., and Jaspers, E. J. ], Universa, Wettere. [35] Wetzel op. c i t . 353 [36] i b i d . 447,448 [37] i b i d . 485 [38] i b i d . 478 [39] Lanczos, C., Applied Analysis, Prentice Hall (1956), 453 [40] Ahlberg, J.H., Nilson, E.N., and Walsh, J.N., The theory of splines and t h e i r a p p l i c a t i o n s . Academic Press, (1967) [41] Pletcher, R., FORTRAN subroutines for minimizing by guasi-newton methods. Research group report. Theoretical physics d i v i s i o n . Atomic Energy Research Establishment, Harwell, England. 148 [42] Fletcher, R., A new approach to variable metric algorithms, Comp. J. 13:317-322, (1970) [43] Powell, M.J.D., and Fletcher R., On the updating of LDL ' f a c t o r i z a t i o n , AERE Harwell Report TP 519 (1972) [44] Fletcher, R., and Powell M.J.D., A r a p i d l y convergent descent method for minimization, Comp. J. 6:163-168 (1963) [45] Gear, B.C., Numerical i n i t i a l value problems i n ordinary d i f f e r e n t i a l eguations. Prentice H a l l , (1971) [46] S t o l l , R.R., Linear algebra and matrix theory, Dover, New York, (1952) 258 [47] Domsellaar, B. van, and Hemker, P.W., Nonlinear parameter estimation i n i n i t i a l value problems. Mathematical Center Report NW18/75, Amsterdam (1975) [48] Federenko, A.Y. And Swift, M.C., Comparative bioloqy of Chaoborus americanus and Chaoborus t r i v i t t a t u s i n Eunice Lake, B r i t i s h Columbia, ' Limnology and Oceanography 17:721-730, (1974) [49] Federenko, A.Y., Instar and species s p e c i f i c diets in two species of Chaoborous, Limnology and Oceanography 20:238-249, (1974) [50] Wetzel op. c i t . 482-483 [51] McKeown, J. I . , Specialized versus general-purpose algorithms for minimizing functions that are sums of squared terms. Math. Prog. 9:57-68 (1975) [52] S u t c l i f f e , V.H.jr., Drinkwater, K., and Muir, B.S., Correlations of f i s h catch and environmental factors i n the Gulf of Maine, J . Fis h . Res. Board. Can. 34:20-30, (1977) [53] Northcote, T.G., and Walters, C.J., Response of zooplankton communities to f i s h introduction, ?erh. Int. Verein. Limnol. (in press) [54] Wetzel op. c i t . 442 149 APPENDIX I General Bounds On States The model represented by eguations 1.1 to 1.9 has the following p r o p e r i t i e s : 1.1 Non-negative Biomasses We f i r s t wish to be assured that for p o s i t i v e i n i t i a l biomasses, we are spared the absurdity of solutions giving negative biomasses. Any tr a j e c t o r y leaving the p o s i t i v e region must cross at l e a s t one of the planes N = 0 x ± = 0 . 1 , 1 But on a l l of these planes, the derivative of the state which has reached zero i s i d e n t i c a l l y zero from eguations 1.1 and 1.2 therefore no t r a j e c t o r y can cross the boundaries of the p o s i t i v e region and a l l biomasses remain non-negative. 1.2 Upper Bounded Nutrient Since x±(t) > 0 and N (t) > 0 from i . 1 , and a± (v) > 0 from eguation 1.7, the rate of change of nutrient dN/dt = N(r-bN - £ a x ±) X * 2 s a t i s f i e s the i n e g u a l i t y dN/dt <_ N(r-bN) = Nb(r/b-N) i " 3 The maximum N= max r(v)/b(v) e x i s t s since r(v)/b(v) < rmax/bmin v > 0 from i n e g u a l i t i e s 1.5 and 1.6. Therefore 150 dN/dt < Nb(N - N) i.4 Integrating, ire obtain i.5 t Where 0 < b m i n t < / b (t) dt < b ^ t and for N (0) > N o i.6 N(t) < N(0) And for N (0) < N i.7 N(t) <_ N Therefore N (t) < max{N(0), N) Since we s h a l l be normally concerned with steady state solutions f o r which N(0) < N, the maximum N ever sustained, then the bound St N = max r(v)/b(v) i.9 w i l l be applicable i n these cases. i.3 Upper Bounded Total Biomass The t o t a l biomass, assuming consistent biomass units i s ° i-9 B = N + I x i The rate of change of t o t a l biomass i s which sub s t i t u t i n g for the biomass rates 1.1 and 1.2 yields n dB/dt - N + Y x. i.10 i-11 dB/dt = N(r - bN - J a i + ^ X i ( e l a i N " m i * 151 How consider the components of dB/dt. The production, P of biomass takes place i n terms P = N(r - bN) i.12 Extrema of production are attained where dP/dN = 0 = r - 2bN i . 13 that i s when N=r/2b which i s a maximum since = -2b < 0 dN2 Therefore maximum production i s given by P where ~> i max 2 P = max { N(r(v) - b(v) N) I N = } = max jr44- i " 1 5 max v 2b (v) v 4b(v) Loss of biomass occurs through r e s p i r a t i o n of consumers which by ineguality 1.9 s a t i f i e s i . 16 - y m. x. < -m . y x. ** i i — min u l Conversion losses , which through i n e g u a l i t y 1.8, s a t i s f y i.17 -N I a. x. + N I x ± e i a. = N I x ± ^ ( e . - l ) < 0 I n e q u a l i t i e s i.16,i.17,i.18 and equation i.10 therefore imply dB/dt < P - m , I x. i ' 1 8 — max min u x Therefore t o t a l biomass cannot increase i f i.19 7 x. > P /m . u i — max min 152 Let B = max(£ x (v), P /m ) + N i-20 u i max mm where N i s the upper bound on nutrient biomass. This bound i s s a t i s f i e d i n i t i a l l y , since i 21 B(0) = I x±(0) + N(0) <_ B For B to exceed B, dB/dt must be posit i v e at B(t) =B for some time t . But i f B(t) =B then from 1.20 T x. + N = max {I x.(v), P /m . } + N i-22 L i I max mm and since N < N subtracting N from both sides y i e l d s J x. > max{y x.(0), P /m . } i.23 L I — ^ l max mm and t h i s implies I x > P /m . i.24 u 1 — max mm And by eguation 1.19 we have dB/dt < 0. Therefore B cannot exceed B as given i n eguation i.20. As i n the bound on nutrient, we w i l l be generally considering steady state so that we may use the expression f o r B which does not involve the i n i t i a l conditions for the biomasses and the steady state bound B on t o t a l biomass i s i.25 Sv B(t) < B = P /m . + N — max min A looser bound i n terms of the environmental response function i n e g u a l i t i t e s i s , since P m a x < T ^ / ^ ^ , r r i ?6 „/ \ max r max . -, \ min mxn 153 i.n Positive Biomasses For F i n i t e T Since nutrient has an upper bound N and consumer biomass has an upper bound B, and i n e q u a l i t i e s 1.5,1.6 and 1.7 hold, i.27 SW dN/dt > N(r . -b N - a B) — mm max max Therefore N(t) > N(0) exp {(r . - b N - a B) t} i - 2 8 v ' — . ' r mm max max which i s posit i v e f o r t f i n i t e . S i m i l a r l y , since N(t) > 0, and equation 1.7,1.8 ,and 1.9 hold d x./dt > x,(-m. ) x.29 i — 1 imax Therefore i.30 x(t) _> x(0) exp (-mimax t) Thus biomass cannot become zero i n f i n i t e t . 154 APPENDIX II Deterministic I n s t a b i l i t y Of Hulticonsumer Community He consider a fixed v f o r which the set of nutrient unique species k for which the maintenance requirement i s minimal. I t i s assumed that there w i l l only be a countable number of points i n the ranqe of v where two nutrient maintenance l e v e l s w i l l be i d e n t i c a l , and that these points are i n f i n i t e l y improbable r e l a t i v e to the continuum of possible choices of v and the minimum of i i . 1 i s therefore unique except f o r these values of v. For v constant, functions of v, e i (v) ,a.± (v), m± (v) are constant, therefore the solutions to equations d x./dt = x. (e. a. N - m.) il« 2 i i l l 1 qiven by i n t e q r a t i o n as maintenance lev e l s of consumers, N^v) =mi (v) /e ±(v) a ±(v) has a N = min N. (v) i i . 1 i t i i . 3 x.(t) = x.(0) exp (/ e. a. N - m. dt) i l r . ; o i i i s i m p l i f y to t i i . 4 o Defininq i i . 5 o He have 155 x ± ( t ) = x^O) exp {(e i a ± <N>t - ra±) t} ii.6 The upper bound on biomass requires that x ^ t ) < B for a l l t , which with equation i i . 6 implies that i i . 7 (e ± a ± <N>t - m ±)t ± ln(B/x±(0)) for a l l t > 0 Dividing by t and takinq the l i m i t on both sides as t ~* °° the ine q u a l i t y e, a. <N> - m. < 0 8 i i 00 l — holds for a l l species at steady s t a t e . This must be true for species k, but from i i . 1 i i . 9 V ek \ < m i / e i v 1 * k Therefore a l l species except species k s a t i s f y the s t r i c t i n e q uality i i . 1 0 eJ a, <N> - m. < 0 i ^ k i i 00 i A persistent solution x(t) i s defined as x ±(t) > e > 0 for a l l t. i i . 1 1 Therefore a persistent species must s a t i s f y the condition An (—ITVT) < (e J a, <N> - m.) t for a l l t > 0 i i . 12 x^O) — i i t i ' Takinq the l i m i t as t -*• 0 0 as before, the i n e q u a l i t y 0 - e i a i < N > » " m i i i . 1 3 must be s a t i s f i e d for species i to p e r s i s t . But only species k 156 can s a t i s f y the persistence condition. Therefore the species k for which Nfc(v) = min mi(v)/ei(v) a^v) i i . 1 4 i s the only species which can s a t i s f y the necessary condition fo r persistence i i . 1 3 and the bounding condition i i . 8 . For species k to pe r s i s t , i t i s necessary from combining i n e q u a l i t i e s i i . 8 and i i . 1 3 that e^a^ <N> -mk=0. This condition can be s a t i s f i e d i f and only i f r/b > \ / e ] / i \ r since species k cannot survive i f i t s requirements are not below the maximum attainable nutrient l e v e l r/b. Solvinq for the equilibrium of the one consumer, one nutrient system, we need only consider non-neqative equilibrium values. Two solutions are possible from posi t i v e i n i t i a l conditions. Solvinq f o r the stationary points of equations 1.1 and 1.2, i i . 1 5 N = r/b x = 0 and N = V V ^ k x = (r - b r n ^ a k)/a k i i - 1 6 The second solu t i o n i s obtained when r/b > mi/e^a^, the f i r s t when t h i s condition i s not s a t i s f i e d . A l i n e a r s t a b i l i t y a nalysis of equilibrium 2 w i l l show i t to be asymptotically stable i n the small when x >0, and since i t i s the only s i n q u l a r i t y i n the reqion bounded by separatrices N=0,x=0, a l l t r a j e c t o r i e s converge to t h i s equilibrium i f i t i s in the 157 p o s i t i v e region. The f i r s t s olution i s stable i f the second soluti o n i s not f e a s i b l e .
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Coexistence of species in a fluctuating environment
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Coexistence of species in a fluctuating environment Fitzpatrick, Gordon James 1977
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Coexistence of species in a fluctuating environment |
Creator |
Fitzpatrick, Gordon James |
Publisher | University of British Columbia |
Date Issued | 1977 |
Description | A dynamic model in which multiple consumers of a single nutrient may coexist in a fluctuating environment is given. Only one consumer can persist in a fixed environment, but coexistence may be produced by effects of a fluctuating environmental variable on nutrient utilization differing between consumers. An approximate solution is given for the non-autonomous Lotka-Volterra-Verhulst ordinary differential equations of the model together with heuristic sufficient conditions for construction of a persistent multispecies consumer community. Computational examples demonstrate persistence of an idealized example community for periodic and random environmental fluctuation. Two further examples demonstrate that environmental fluctuation can produce coexistence when environmental variables, standing crops, assimilation efficiencies, primary productivity, utilization rates, and respiration rates are comparable to a tropical grassland, and an oligotrophic temperate lake. The sensitivity of model solutions to functional variations of the component species may be rapidly and accurately calculated. This allows the identification and estimation of unknown species functional responses from time series data of biomasses and a measured environmental variable. Unknown functions of an environmental variable are approximated by a Tchebycheff polynomial expansion in that variable. Unknown coefficients of these expansions are the parameters of the model. These parameters are determined by the unconstrained minimization of the squared deviations of the logarithm of biomass observations and model differential equation solution using a Quasi-Newton algorithm. This least squares estimator was applied to a one year biomass time series of four zooplankton grazers, phytoplankton, and average lake temperature of a small oligotrophic lake. Application of the model to this grazer zooplankton community gives evidence of partial stabilization due to environmental fluctuation in a natural community. It is concluded that environmental variation, which is often assumed on theoretical grounds to be destabilizing, should rather be considered as one of the bases of community persistence. |
Subject |
Variation (Biology) Animals - Variation |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-02-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0065003 |
URI | http://hdl.handle.net/2429/20671 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1977_A1 F58.pdf [ 6.87MB ]
- Metadata
- JSON: 831-1.0065003.json
- JSON-LD: 831-1.0065003-ld.json
- RDF/XML (Pretty): 831-1.0065003-rdf.xml
- RDF/JSON: 831-1.0065003-rdf.json
- Turtle: 831-1.0065003-turtle.txt
- N-Triples: 831-1.0065003-rdf-ntriples.txt
- Original Record: 831-1.0065003-source.json
- Full Text
- 831-1.0065003-fulltext.txt
- Citation
- 831-1.0065003.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0065003/manifest