A Methodology for Mathematical Modeling of Biological Production -by R. P. Canale Department of Civil Engineering The University of Michigan Report.to The University of Mich i gan Sea Grant Project Covering. Jan. 1970 Jan. 1971

(em kS <J., w An. IN e*R ^.X '',>... 0< ~\o

INTRODUCTION In its Lake Erie Report, the FWPCA concluded that the principal reasons why Lake Erie has a pollution problem are the neglect of citizens and governments, and the lack of knowledge concerning the dynamics of lake ecology and the causes of lake pollution. An expert knowledge of an ecosystem can only be developed through the extended work of professionals from a variety of disciplines. Water resource planners, managers, and politicians, who are not ecosystem experts, must make decisions involving large amounts of public money regarding the impact of waste disposal. Obviously, a workable model capable of simulating natural systems and quantitatively predicting the natural population dynamics resulting from various degrees of pollution would be valuable in water resource planning. These ideas were also expressed in the original Sea Grant proposal submitted to NSF in 1969. Figure I is taken from this document, and illustrates the roles of modeling in water resource management. Figure 2, also from this document, shows the various sub-models which must be used in the synthesis of a viable decision making tool. This report discusses the results of studies concerned with developing the methodology required for Implementing biological production simulations using mathematical methods. Several workers have presented arguments for the utility of the mathematical modeling approach in ecology. Particularly noteworthy are the publications of Rapoport (1966), von Bertalanffy (1962), Waterman and Morowitz (1965), and Watt (1962). As an example, we know that the activities of our increasing population have accelerated the aging of the Great Lakes at a time when more human beings than ever are dependent on them. However, models are not available to balance the conflicting uses and thus managing the resources of the

FUNCTIONS OF THE INTEGRATIVE CONSORTIUM Ill Definition of Environmental Tolerance Limits IV. Definition of New Concepts Resource Systems 11 / / nManagement objectives VI....i~~~~~~ -, - 1 Decision Systems m 1~Im I I/| ~ Decision Systems |~ II*^ [ ~ Application of Current Data Aquisition L Systm Knowledge & Predictions Input to regulation gt& qIIto Steto Do of uses and managment Assembly M odcllking V Decisionmaking of attributes Optimization on Basis of!') \ \ Multiple Costs Benifits VI Evaluation of Future Social & Economic Impacts VII Evaluation of Future Technological Impacts Figure I —Relationships among the Integrative Functions of The University of Michigan Sea Grant Program.

Cn w CL z 'SOIL C1) DISSOLVED SOLIDS SUSPENDED SOLIDS MICROFLORA z Z O 1: MACROFLORA L) DISSOLVED SOLBOTTOM SEDIMENTS / HERBVOR SH I CURRENTS ATMOSPHERE TRIBUTARY WATER a- RRENTS ~ STREAMS WITHDRAWALS RECOVERY RESOURCES H DISSOLVED SOLIDS GASES MIGRATORY MUNICIPAL SEDIMENTS FISH a- SUSPENDED SOLIDS FAUNA INDUSTRIAL BEDROCK INVERTEBRATES H BIOTA AGRICULTURAL OIL & GAS BIRDS O MAMMALS Figure 2 — Components of the Subsystem of Marine Resources of Traverse Bay, with Input and Output Channels. 3

Fish,A~gs.. rotozoa Bcteria Inorganic &lOrganic Material Figure 3. Simplified food web. 4

Great Lakes system is difficult. We hope that the work described herein is a step toward the development of models useful in such decision-making processes. The microbial community in a lake is affected by the amount and type of organic and inorganic nutrients in the aquatic environment. The structure of a typical community is fairly well-defined and results from the flow of energy through a diverse network called a food web. The food web is illustrated in Figure 3. It is this consistency of structure that enables a vast natural system to be stimulated by a much less complicated model. The food web concept provides a framework for a rational interaction between ecologists and engineers concerned with population dynamics of large, complex, and highly variable ecosystems and water resource planners concerned with the impact of waste effluents on our lakes. The synthesis of a model of a complex natural ecosystem requires the development of mathematical descriptions of the diverse chemical-biological interactions that occur. Researchers have devoted considerable efforts to obtain mathematical formulations for the primary production of photosynthetic green plants and subsequent inorganic nutrient consumption. Others have modeled the decomposition of dissolved organic material by the heterotrophic microorganisms. Quantitative mathematical models for describing the predation effect on algae and bacteria by the protozoa, the lower invertebrates, and insects have received comparatively little attention. This is true despite the fact that these predators are intimately associated with the overall eutrophication picture. The predator food chain has the effect of transferring, through a series of steps, the energy contained in the nutrient substances to the successively higher forms of aquatic life with a loss of potential energy occurring at each transfer step. Thus, the 5

amount and kinds of aquatic life in a lake are determined by the complex interrelationship of water temperature, sunlight intensity, and nutrient concentrations on one hand, and the grazing action of various predators on the other. The report by the Subcommittee on Environmental Improvement, Communittee on Chemistry and Public Affairs, American Chemical Society has recommended that further research should be directed to developing improved mathematical descriptions of predation in such natural food chains. METHODOLOGY Model Synthesis Mathematical models could be used to aid in a rational approach to the management of diverse water resource systems. The need for a completely rational procedure becomes more acute as our utilization of the environment increases. Mathematical models of biological production can be developed as proposed in Figure 4. Kinetic or rate expressions for the growth of each biological species are required along with the associated nutrient utilization rates. The specific growth and nutrient uptake rates are dependent on the concentration of nutrients, temperature, pH, and light intensity. Additional factors include various microbial interactions such as predation, competition, and symbiosis. Information of this type must be expressed in mathematical form to be of maximum usefulness. Most natural bodies of water have complicated and unique internal flow patterns and mixing characteristics. These mechanisms interact with a diverse network of productivity zones and stratifications. Together these factors affect 6

Biological Kinetics Flow Patterns, Material and Rate Processes Mixing, and Balance Circulation j Synthesis of _ Field Plankton Model Observations....1~^~ -.....Performance Design of Predictions Laboratory Experiments Comparison of Field Data and Model Predictions Laboratory Results Figure 4. Dynamic modeling procedure. 7

the rates of biological production and nutrient uptake. It is proposed that a given body of water be divided into several zones of relatively uniform composition. The resulting uniform cells then represent a three-dimensional network which can be used to account for various phenomena. For example, horizontal sectioning may divide a lake to account for shoreline behavior, bays, and circulation zones. Figure 5 illustrates the procedure for a case like Lake Michigan. Vertical sectioning can account for transfer of gases from the atmosphere, variation of sunlight with depth, temperature and oxygen distributions, sedimentation of detritus, migrations, and nutrient releases from bottom deposits. Figure 6 shows a possible cell structure to account for the various actions. Each relatively uniform cell interacts with adjoining cells and a transfer of material among the cells results from convective transport, dispersive transport, migration and sedimentation. Each of the uniform cells constitutes a unit which can be conveniently described by writing a non-steady material balance. Within each cell the behavior of any chemical or biological species i can be described as,, Rate of Change of i Rate of Input of i Rate of Output of with time within a by convection, dis- i by convection, cell persion, sedimenta - dispersion, sed(weight/time) tion, or migration imentation, or from adjacent cells migration from (weight/time) adjacent cells (weight/time) Rate of Production Rate of Disapof i by growth, pearance of i by -+ excretion, or - uptake, predation, dissolution within respiration, death, cell or precipitation (weight/time) (weight/time) (I) 8

Zone 1 Zone 2 Figure 5. Hypothetical horizontal sectioning of Lake Michigan into uniform zones. 9.

Sunlight Exchange of Exchange of Oxygen Carbon Dioxide Zone I V I Zone 1 Zone 2 Sedimentation M igrations \^ ~Zone \ /y^ ^^^ ~~~~Sediaentatlon ^^ ^^^^ l ~~~~~~Oxygen ^ I^Zone - ^ ^ ^ Organic Material Benthos Figure, 6. Vertical sectioning of a lake into uniform cells.

This method of accounting for the spatial distribution of species, nutrients, temperature, and light Is useful In a large variety of problems, while at the same time having the advantage of mathematical simplicity. The resulting differential equations are ordinary rather than partial and therefore are an alternative to the earlier models of Riley, Stommell and Bumpus (1949), Kierstead and Slobdkin (1953), and Strickland (1960). The reported work is a generalization of the approach of Steele (1958) where a thermally stratified sea was divided into two homogeneous layers. The proposed finite element approach is extended to include horizontal as well as vertical transport. Patten (1965) has utilized a similar idea for elemental volumes. As a result, his equations were partial rather the ordinary and require a continuous integration over time and position. Obvious difficulties arise when dealing with irregular shapes or unusual boundary conditions. The immediate benefit of having a mathematical model incorporating kinetics, flow patterns, and mixing is the availability of a theoretical basis for designing pilot or feasibility experiments. Thus it is possible to facilitate the design of experimental equipment, the determination of the most significant variables to be measured in an experimental program and the interpretation of experimental results. Laboratory results and field data can be used via the mathematical model to construct a full scale model of an aquatic ecosystem. The full scale model then permits a prediction of the system behavior which can be compared with survey data. A comparison of the prediction and the survey data can result in feed-back to the synthesis step resulting in model alteration. Thus, the modeling procedure is viewed as a II

dynamic process of synthesis, prediction, comparison, and modification. A Test Ecosystem A useful predictive model of the ecological processes involved in aquatic systems must clearly be limited in scope. It is not expected, required or even desirable in many cases to attempt to describe all the possible interactions. However, a useful model must contain enough detail to predict the gross behavior or characteristics of the system in response to outside influences or driving forces. In order to test the feasibility of the reported techniques for modeling, a simplified model of a uniform cell has been selected as shown in Figure 7. The major constituents of the biology of the cell are Illustrated in Figure 8. Inorganic nutrients such as phosphorus and nitrogen are utilized by heterotrophic bacteria parallel with the oxidation of organic carbon. These same Inorganic nutrients are required by the algae with the utilization of inorganic carbon as stimulated by energy from sunlight. The interaction of algae and bacteria is one of both competition and symbiosis. The bacteria produced during the destruction of organic carbon will lead to the development of a population of predator-x. The algae will stimulate the growth of a different class of organisms, predator-y. Either predator-x or y will result in predator-p. The photosynthesis and respiration of the community constitute various sources and sinks of dissolved oxygen and carbon dioxide. Overall, the test system neglects many features of a natural population, but does contain enough flexibility to hope that a useful simulation of the population dynamics can result. An obvious factor which the model does not 12

Uniform Volume N X B Y C P /1(10 A B1 '1 1 Figure 7.Uniform A l Figure 7. Uniform rectangular volume.

0 CIA NUTRIENTS C02 ^ PREDATOR-P-02 I42~10 OC Figure 8. Microbial and chemical interactions to be simulated with mathematical models.

include are effects of nutrient release by anaerobic bottom deposits. Thus, the results of this work apply to only the euphotic zone. A) Bacteria Growth Kinetics Numerous experimental and theoretical studies have been conducted on the kinetics and stoichiometry of various subsets of the system depicted in Figure 8. Monod (1942) and Herbert et al. (1956) have investigated bacteria growth and organic carbon removal using continuous culture techniques. These studies have resulted in the use of Michaelis-Menton kinetics for bacteria growth, C. rb = 11 (C)b = K ) b (2) b K + C c where rb is the growth rate in mg/I/day, p is the specific growth rate, p is the maximum specific growth rate when C -+, C is the concentration of organic carbon, b is the concentration of bacteria, and K is a saturation constant equal to the concentration of organic carbon when rb/b is equal to "/2. Associated with the bacteria growth is an organic carbon substrate removal whose stoichiometry is defined by the Equation (3), r = - r (3) c Y 'b where r is the rate of carbon removal in mg/I/day and Y is known c as the yield constant and is the ratio of the weight of bacteria produced to the weight of carbon consumed. Equations (1), (2), and (3) can be applied to a homogeneous and well-mixed cell to give, 15

db -^ Q dt + Cb - b (4) dt1 K +0 V.d 0 Q C C dt V CI C) Y 'K + C b () ' C where CI is the concentration of organic carbon in an incoming stream of flow rate Q Into a volume V as illustrated in Figure 9. Bacteria growth in continuous culture with multiple nutritional deficiencies has been considered by Cooney and Mateles (1969). The growth of Aerobacter aerogenes was characterized as a function of the concentration of organic carbon, nitrogen, and phosphate. Batch cultures of a mixed population of bacterta. have been studied under multiple nutrient limitations by Schaezler (1969). The modeling reported herein will attempt to account for bacteria growth with multiple nutrient restriction by introducing a modification of the Monod growth model, C n r = P (Cn) b + + n b (6) c nb where n is the concentration of a restricting nutrient such as phosphate and Knb is a saturation constant. The variation of rb with temperature can be represented by the Arrhenius function. B) Algae Growth Kinetics There have been numerous studies dealing with mathematical models of microscopic photosynthetic plant growth. The major variables affecting the growth are nutrient concentration, light intensity, and temperature. However approaches vary: Patten (1968) 16

Mixing Concentration Flow Rate Q / Bacteria Produced = b Concenti o f Remaining/ Organic Concent^^io^ F Remaining Organic Carbon = C Carbon c Figure 9. Continuous culture of bacteria utilizing organic carbon.

has published, an extensive review of plankton production models and in each case where nutrient concentration was considered, it was assumed that r c n (7) a where r is the growth rate of algae in mg/l/day. Others, such a as the Joint Industry-Government Task Force on Eutrophication (1969), have suggested the use of Michaelis-Menton kinetics to express the growth rate of algae as a function of nutrient concentration, r = (n) a = (K n a a K ' n na where ~ Is the maximum value of the specific growth rate p, a is the concentration of algae, and K is a saturation constant. na Dugdale (1967) presented a theoretical framework for nutrient limitation in the sea and also used Michaelis-Menton kinetics for plankton growth. Several authors have discussed the variation of algae growth with light intensity as influenced by water depth and self-shading. In general, the modeling of such effects is complicated, and the mathematics becomes very difficult, e.g., see Sverdrup (1953), Tailing (1957), Strickland (1960), Murphy (1962), Steele (1962), Patten (1965), and Vollenweider (1965). The homogeneous finite element cell approach, as suggested in this report facilitates the mathematical analysis as simply as possible. Within a uniform volume, the light intensity is a function of time, as Illustrated In Figure 10. Photosynthesis can be 18

Average Intensity of Sunlight -itthin Uniform Volume I max Time 6 a. m, 7 p.m. 6 a.m. Figure 10. Hypothetical variation within a uniform volume as a function of the time of day.

related to intensity by a hyperbolic function, as suggested by Vollenweider (1965). There is an inhibition interaction between the plankton population and photosynthesis. Probably the most important effect is the self shading of the algae and the shading due to the bacterial population. This report explores the use of an exponential Inhibition function to describe this shading decrease of available light energy, Available light = I(t) exp[-k(b+Oa)] (9) where k and B are constants. Combining Equations (8) and (9) and introducing a saturation relationship for light results in a model for growth due to photosynthesis, ra " )(K( n I(t)) -k(b + Sa) ) r K + K e+ln))e a (10) a K+V K ~+ it e na I where K is a saturation constant for light. The kinetic expression given by Equation (10) is sufficiently general to account for many phenomena observed during algae growth. Notice that high light intensity inhibition can be included in the model in any finite element by calculating a lowered average ~ reflecting this inhibition. The variation of ra with temperature can be described with Arrhenius function. The kinetic expression given by Equation (10) is considerably more realistic than most previously published models. For example, Parker (1968) assumed a linear variation between r and nutrient concentration and light intensity and a did not account for shading. The assumption of Parker can result 20

in serious misinterpretations of experimental data and can result in unrealistic model behavior. The consequences of different kinetic assumptions in microbial interaction modeling has been discussed by Canale (1970). C) Predator Growth Kinetics The classical work of predator growth kinetics was carried out in 1925 and 1926 by Lotka and Volterra. These investigators studied models where predator specific growth rates had a linear variation with the prey concentration. Researchers such as Davidson and Clymer (1966) and Parker (1968), to name a few, have retained these historical assumptions in plankton modeling. Others such as, Canale (1968), Bungay and Bungay (1968), Drake et al., (1966), and Curds and Cockburn (1968) have used MichaelisMenton kinetics for predator growth. Canale C.1970) has discussed the importance of careful consideration of kinetic expressions in predator-prey models. This report investigates the useful ness of Michaelis-Menton kinetics for predator growth as expressed by Equations (II), (12), and (13) r m(b)x = m ( b ) x (II) x Kb +b b a r = 6(a)y =6 (K a+) y (12) y K a a ^ x + ay ) (13) r = y(x,y)p = y (K P ~p ~ ~x + x + ay 2 xy 21

where rx, ry, and r are the growth rates of predator x, y, and p and where m, 6, y are the maximum values of the specific growth rates m, 6, and y for predators x, y and z and Kb, Ka, K and a are constants. The abili ty of predator p to utilize either x or y is described according to Equation (13), and has not been tested experimentally. The temperature affect on the predator growth can again be included by the Arrhenius equation. D) Simulation of Example Problem The general material balance equation (Equation I) and the kinetic expressions given by Equations (6), (10), (II), (12), and (13) have been utilized by Canale (1970b) in a preliminary project to construct a plankton model. The results are given by Equations (14) - (22), dn Q ra rb dt V I - n) - Wa - Wb (14) da Q (14) r dt V (al - a) + r a - -y (15) dt (Y - Y) + Wp x + ay db_ x (b (- b b) + r(17) dt V b WX dx Q r__x (.18) dt (x - x) + rx - -- ( ay d- C(C - C) WC (19) dt V I - Wb rb 22

t V (p -p) + r - (mining) (20) dC02 Q WC02 Production by respiration, dt V (C02 - C02) W ra + atmospheric absorption and dt V I Wa a chemical equi libria (21) d02 = Q (02 - 02) - community respiration + photosynthesis + atmospheric reaeration (22) where, the subscript I refers to the input concentration of the various species from adjacent cells, and Wa is the ratio of the weight of algae produced per weight of nutrient utilized, Wb is the ratio of the weight of bacteria produced per weight of nutrient utilized, Wy is the ratio of the weight of y produced per weight of algae utilized, Wp is the ratio of the weight of p produced per weight of y or x used, WC02 is the ratio of the weight of carbon dioxide required.per weight-of nutrient utilized, Wx is the ratio of the weight of x produced per weight of bacteria utilized, and Wc is the ratio of the weight of organic carbon required per weight of nutrient utilized. In order to keep the analysis relatively simple, the dispersion, sedimentation, excretion, and respiration terms have been considered negligible, and constant temperature has been assumed. Yet, because of the complicated form of the various growth and removal rates, the mathematics are still complex. The solution of the above system of nonlinear ordinary differential equations requires the use of numerical or analog methods. The transient behavior of the test model has been studied with analog techniques by Canale (1970a). In that study typical values of 23

the various parameters were obtained from published reports, or more typically estimated from intuition. Table I shows the estimated parameter values used in the computer program. Figure II illustrates the results of the computer output for a typical run. The high initial nutrient concentrations are quickly removed by the bacteria and the algae. The bacteria and algae result in the development of significant numbers of predator x and predator y. Due to the coupling of the equations and resulting feed back the growth of predators-x and y cause the nutrients to be removed at a reduced rate as evident from the change in the shape of the nutrient output curve. The predator-p gradually increases at the expense of x and y. This, by way of the feed back loops, begins to have a noticeable effect on the bacteria and organic carbon concentration. The plausible results of these preliminary experiments are encouraging, and warrent further investigation as suggested in this report. Furthermore, the experimental results of Canale (1968) and Curds and Cockburn (1968) for an organic carbon-bacteria-protozoa subset system are additional impetus justifying further study of the reported methodology. The properties of the solution of systems of nonlinear differential equations depend on the magnitude of the eigenvalues as reflected by the values of the parameters. It was not surprising therefore, to observe damped and steady state oscillations in the system as the parameter values were changed in several computer runs, It is anticipated that digital-analog hybrid computers may be 24

SATURATION CONSTANTS MAXIMUM SPECIFIC GROWTH RATES Kc = 5 mg/l = 3 day Kb = 20 mg/l M = 1 day 1 K =.5 mg/l y =.1 day na I.1 Imax Knb =.5 mg/l 6 = 1 day 1 Ka = 15 mg/l Ka =15 mg/ FLOW RATE AND VOLUME Kxy = 5 mg/l 3 xy =5- 5mg/i -2Q = 2 x 107 ft3/day 9 3 STOICHIOMETRIC COEFFICIENTS V = 10 ft W 50 C02 INPUT CONCENTRATION Wa = 30 CI = 50 mg/l Wy =.5 Wp =.5 aI 3 mg/l Wc = 100 YI = 1.5 mg/l Wb = 60 bI = 6 mg/l Wx =.5 xI = 3 mg/l PI =.75 mg/l TABLE 1. ESTIMATED PARAMETER VALUES FOR COMPUTER SIMULATION 25

)1,, I~~ I lT q~~~~~~~~~~~~~~~~~ r0 ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~r i / ~~~~~~~~~~~~ (L~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~) o O 40 *~_~_ ~ ~ CON>TR AT ION OF /E C IES ^^^~~^*^~~~ f \ ^" lJ ~ CONOENTRATION~~~~~~~~~~~~ ~ OFSECE / / \ ^ eg ^~~~~~~~~~~~~~~~~~~P ^^ / \ ^ ^ ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~tip, /\ ' ca~~~~~~~~~~~~~~~o ^^ ^ / \ ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~,z ~ \h; Q t sl ^K \ Q>~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ^^ ^ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~) - ^ s. \X:$ '^ (^ ^ I, \ <^r \ ^ 1^^~~~~~~~~~~~~~4 CONCENTRATION OF SPECIES

used to advantage as illustrated in Figure 13. RECOMMENDATIONS The specific objectives of our work are concerned with developing ecosystem modeling techniques using as a starting point the model given by Equations (14) to (22) and moving toward more general ecosystem models for water resources planning and management. Laboratory culturing, mathematical analysis, and hybrid computer simulation are proposed to explore its potential. Experimental Studies The experimental work suggested herein consists of a network of continuous culturing studies. The cultures will start simply and develop gradually and become more complicated. The first efforts should concentrate on continuous culturing of a mixed population of bacteria serving as food for a mixed protozoan population, as illustrated by the food chain, Organic carbon __q heterotrophic __. holozoic bacteria protozoa The experiimental procedures and equipment for such culturing have been described by Canale (1968) and are currently being used in his laboratory. Parallel to these efforts, it is suggested that the symbiotic interrelationship and competition between mixed algae and bacteria populations be studied using continuous culturing techniques. The potential importance of these interactions has recenlty been pointed out by Kuentzel (1969), and Kerr et al., (1970). These symbiotic interactions are illustrated in Figure 12. 27

ORGANIC CARBON OXYGEN CARBON ~ ALGAE BACTERIA ~ DIOXIDE ~ DIOXIDE NUTRIENTS Figure 12. Symbiotic interaction between bacteria and algae. 28

uJ LU Q'ATA STORE TRANSFER ANALOG UPDATED RAW INFORMATION COMPUTER ARAMETER TO PARAMETER CONTROL VALUES VALUES DIGITAL COMPUTER TRANSIENT SOLUTIONS MODEL PARAMETER FROM EQUATIONS (14D-(22) VALUES Figure 13. Hybrid computer system for ecological simulation.

Additional experimental work should concentrate on extending the continuous culture work to more complicated systems. Specifically, the predation of Daphnia on the algae population used in the symbiosis studies should be investigated in continuous culture. The overall response could be simulated with Michealis-Menton kinetics although in this case the characteristics of the Daphina are relatively more complicated. Of course, studies have been conducted by many investigators on Daphina nutrition and growth and these results may be useful. The reports of Hall (1964), McMahon and Rigler (1964), Smith (1963), and Wright (1965) should prove particularly useful. Here the need for a dynamic approach to modeling with continuous modification is acute. An ecological model with the degree of sophistication suggested here must be thoroughly tested under controlled laboratory conditions before there is any hope of predicting ecological events in a natural environment. After such laboratory studies have been completed it should be possible to analyze existing field data and propose future in situ experimental programs with a reasonable chance to achieve the desired objectives. Analytical and Numerical Studies Models for plankton production almost invariably result in systems of nonlinear differential equations. A notable expection to this pattern is the linear model of Patten (1966). The mathematical theory for systems of nonlinear differential equations is not as fully developed as for systems of linear equations. Existence and uniqueness theorems only apply to a limited class of cases. Furthermore, 30

the behavior of the solution of nonlinear equations may be quite complicated and phenomena can exist which cannot possibly result from linear systems. As an example, only nonlinear systems can result in "limit cycle" periodic solutions. The behavior of the solution of nonlinear predator-prey models has been studied by Canale (1970) using classical methods of nonlinear mechanics. The techniques of phase-plane and singular-point analysis proved particularly useful. It is worthwhile to analyze a nonlinear model such as Equations (14) to (22), or a model with experimental modifications, using the methods of nonlinear mechanics. The singular points of the model can be determined using the suggestions of Stern (1965). The eigenvalues near each singularity can be expressed as a function of the parameters. This will permit the characterization of the solution near the singular points. Thus, after the system parameters have been evaluated experimentally, it should be possible to state qualitatively the expected variation in the dependent variables as a function of time without the need to develope solutions numerically. The resulting savings in computer time will be substantial in large systems. Furthermore, the iegenvalues give direct information concerning the stability of the system. Stability information is of course vital before rational decisions can be made concerning resource allocation and distribution. It is also suggested that the useful ness of a hybrid computer system be explored as a tool for simulating ecology phenomena. Our analog computer is an Applied Dynamics-4 which has 100 amplifiers and 24 multipliers. The digital is a PDP-9 and has a core storage of 18K. The analog is an ideal integrator, and can be used to obtain rapid dynamic solutions of Equations (14)-(22). 31

The digital computer has the attractive features of its ability to store and retrieve data, and to perform logical and algebraic operations, Figure 13 shows the expected interaction between the analog and the digital computer in the hybrid system. Field generated raw data on species changes, temperature, circulation, and sunlight are converted to values of model parameters and updated values are sent to the analog for simulation. The solutions of the models as presented in Figure II are for the case of a single homogeneous cell. It is expected the realistic problems would involve the use of at least twenty cells. However, it should be noted that the capabilities of our analog system were completely utilized for just one cell, thus, it is anticipated that an iterative approach must be adapted. These computing difficulties are highly significant and require considerable expertise in methods of numerical integration. 32

REFERENCES Bungay, H.R. and Bungay, M.L., "Microbial Interactions In Continuous Culture," Advances in Applied Microbiology, 10, 269-290 (1968) Canale, R.P., "Predator-Prey Relationships in a Continuously Fed Biological Reactor," Ph.D. Thesis, Syracuse University, Syracuse,' N.Y., (1968) Canale, R.P., "Predator-Prey Relationships in a Model for the Activated Sludge Process," Biotechnology & Bioengineering, 11, 887-907, 1969 Canale, R.P., "An Analysis of Models Describing Predator-Prey Interaction," Biotechnology & Bioengineering, 12, 353-378, 1970 Canale, R.P., Unpublished progress report on productivity modeling for the Sea Grant Program of The University of Michigan, 1970b Cleaning Our Environment, The Chemical Basis For Action, A Report by the Subcommittee on Environmental Improvement, Committee on Chemistry and Public Affairs, American Chemical Society, Washington, D.C., 1969, pages 13 and 102 Cooney, C.L. and Mateles, R.I., "Multiple Nutritional Deficiencies in Continuous Culture," presented at the 158th National American Chemical Society Meeting, New York, September 1969 Curds, C.R. and Cockburn, A., "Studies on the Growth and Feeding of Tetrahymena pyriformis in Axenic and Monoxenic Culture," Journal of General Microbio ogy, 54,343-358 (1968) Davidson, R.S. and Clymer, A.B., "The Desirability and Applicability of Simulating Ecosystems," Ann. N.Y. Acad. Sci., 128, 790-4, 1966 Drake, J.E., Jost, J.L., Fredrickson, A.G. and Tuschiya, H.M., "The Food Chain," Bioregenerative Systems, A conference held in Washington, D.C., November 15-16, 1966 sponsored by NASA and the American Institute of Biological Sciences, NASA SP-165, Library of Congress number 68-60345, U.S. Government Printing Office, Washington, D.C., (1966) Hall, D.J., "An Experimental Approach to the Dynamics of a Natural Population of Daphnia galeata mendotae," Ecology, 45, 94-112, 1964 Kerr, P.C., Paris, D.F. and Brockway, D.L., "The Relationship Between Carbon and Phosphorus in a Mixed Heterotrophic and Autotrophic Population in an Aquatic System," to appear in the Proceedings of the 25th Annual Industrial Waste Conference, Purdue University, 1970 Kierstead, H., and Slobodkin, L.B., "The Size of Water Masses Containing Plankton Blooms," J. Mar. Res., 12, 141-7, 1953 33

Kuentzel, L.E., "Bacteria, Carbon Dioxide, and Algal Blooms," JWPCF, 41, No. 10, 1737-1747, Oct. 1969 Lake Erie Report, A Plan for Water Pollution Control, U.S. Department of the Interior, FWPCA, Great Lakes Region, August 1968, page I Lotka, A.J., "Elements of Physical Biology," The Williams and Wilkins, Co., Baltimore, MD. (1925) McMahon, J.W. and Rigler, F.H., "Feeding Rate of Daphnia magna Strauss in Different Foods Labelled with Radioactive Phosphorus," Limnol and Oceanogr., 10, 105-113, 1964 Murphy, G.I., "Effect of Mixing Depth on the Productivity of Freshwater Impoundments," Trans. Am. Fish Soc., 91, 69-76, 1962 Parker, R.A., "Simulation of an Aquatic Ecosystem," Biometrics, Vol. 24 (4), 803-821, 1968 Patten, B.C., "The Biocoenetic Process in an Estuarine Phytoplankton Community," Oak Ridge National Laboratory, Tenn. (ORNL-3946) 1965 Patten, B.C., "Community Organization and Energy Relationships in Plankton," Oak Ridge National Laboratory, Tenn. (ORNL-3634) 1965 Patten, B.C., "Mathematical Models of Plankton Production," Int. Revue ges. Hydrobiol., 53, 3, 357-408, 1968 Provisional Algal Assay Procedure, Joint Industry Government Task Force on Eutrophication, P.O. Box 3011, Grand Central Station, N.Y., N.Y,, February 1969 Rapoport, A., "Mathematical Aspects of General Systems Analysis," p. 3-11, In L. von Bertalanffy and A. Rapoport (ed.). General Systems, Vol. II, Society for General Systems Research, Ann Arbor, Mich., 1966 Riley, G.A., Stommel, H., and Bumpus, D.F., "Quantitative Ecology of the Plankton of the Western North Atlantic," Bull. Bingham Oceanogr. Col I., 12, 1-169, 1949 Schaezler, D.J., Busch, A.W., and Ward, C.H., "Kinetic and Stoichiometric Limitations of Phosphate in Pure and Mixed Bacterial Cultures," to appear in the Proceedings of the 24th Annual Industrial Waste Conference, Purdue University, 1969 Smith, F.E., "Population Dynamics in Daphnia magna and a New Model for Population Growth," Ecology, 44, 651-622, 1963 34

Steele, J.H., "Plant Production in the Northern North Sea," Pub. Scottish Hom Dept. Mar. Res., 1. 36 pp., 1958 Steele, J.H., "Environmental Control of Photosynthesis in the Sea," Limnol. Oceanogr., 7, 137-50, 1962 Stern, T.E., "Theory of Nonlinear Networks and Systems," Addison-Wesley, Reading, Mass., 1965 Strickland, J.D.H., "Measuring the Production of Marine Phytoplankton," Bull. Fish. Res. Bd. Can., No. 122, 172p, 1960 Sverdrup, H.U., "On Conditions for the Vernal Blooming of Phytoplankton," J. Cons. Explor. Mir., 18, 287-95, 1953 Tailing, J.R., "Photosynthetic Characteristics of some Freshwater Plankton Diatoms in Relation to Underwater Radiation," New Phytol. 56, 29-50, 1957a Vol lenweider, R.A., "Calculation Models of Photosynthesis-depth Curves and some Implications Regarding Day Rate Estimates in Primary Production Estimates," Mem. Inst. Ital. Idrobiol, 18(Suppl.), 425-57, 1965 Volterra, V., "Variazioni e fluttuazioni del numers d' individui in specie animali convienti," Mem. Acad. Lineii Roma, 2, 31 (1926) von Berttalanffy, "General System Theory - A Critical Review," p. I to 20. In L. von Bertalanffy and A. Rapoport (ed.), General Systems, Vol. 7. Society for General Systems Research, Ann Arbor, Mich. 1962 Waterman, T.H. and Morowitz, H.J., (ed.), "Theoretical and Mathematical Biology," Blaisdell Pub. Co., NY.Y., 426 p., 1965 Watt, K.E.F., "Use of Mathematics in Population Ecology," Ann, Rev, Entomol. 7, 243-60, 1962 Wright, J.C., "The Population Dynamics and Production of Daphnia in Canyon Ferry Reservior, Montona," Limnol. and Oceanogr, 10, 583-590, 1965 35

UNIVERSITY UP MIUMIUAN 3 9115 02651 8624111111 3 9015 02651 8624