Global behavior of a generalized aquatic ecosystem model

Ecological Modelling, 45 (1989) 95-110 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

GLOBAL BEHAVIOR OF A GENERALIZED ECOSYSTEM MODEL

9'i

AQUATIC

T. KMEq" i and M. STRASKRABA 2 i Institute of Experimental Biology and Ecolo~', Centre of Biological and Ecological Sciences of the Slovak Academy of Sciences, Bratislava (Czechoslovakia) 2 Biomathematical Laboratory, South Bohemian Biological Centre, Czechoslovak Academy of Sciences, Ceskd Bud~jovice (Czechoslovakia)

(Accepted 7 October 1988)

ABSTRACT Kmei, T. and Stragkraba, M., 1989. Global behavior of a generalized aquatic ecosystem model. Ecol. Modelling, 45: 95-110. This paper investigates steady-state behavior of a simple dynamic model of an aquatic ecosystem (series AQUAMOD, model GIRL). When throughflow is zero, the biological components vanish and only dissolved phosphorus remains. For nonzero throughflow, one asymptotically stable in the large steady state was found for every condition investigated. For periodically varying inputs a periodic trajectory of steady-state solutions for different t was obtained. The results of steady-state investigation are confronted with both dynamic model solutions and natural observations.

INTRODUCTION This p a p e r is d e v o t e d to an e x a m i n a t i o n of qualitative b e h a v i o r of a simplified, generalized aquatic e c o s y s t e m b y m e a n s of s t e a d y - s t a t e solution of its model. A d e m o n s t r a t i o n will b e m a d e of s o m e e l e m e n t a r y p r o p e r t i e s o f a simple e c o s y s t e m u n d e r c o n s t a n t a n d seasonally v a r y i n g c o n d i t i o n s . In addition, an a t t e m p t is m a d e to c o n f r o n t the s t e a d y - s t a t e b e h a v i o r of the m o d e l with the system's dynamics. MODEL DESCRIPTION T h e m o d e l used in this p a p e r is d e r i v e d f r o m the series A Q U A M O D b y Stra~kraba a n d G n a u c k (1985). I n t e n d e d for a G e n e r a l i z e d I m i t a t i o n of Reservoirs and Lakes ( G I R L ) , it consists of three e q u a t i o n s d e s c r i b i n g p h o s p h o r u s as the limiting element, a n d p h y t o p l a n k t o n a n d z o o p l a n k t o n in 0304-3800/89/$03.50

© 1989 Elsevier Science Publishers B.V.

96

a uniformly mixed layer of water with specified depth (ZMIX). It may represent either a shallow lake or the epilimnic upper layer of a deep lake. The reservoir character is due to variable throughflow ( Q / V ) , the additional physical characteristics being represented by the extinction coefficient of water for light (EPS) and annually varying temperature (TEMP). To simulate the annual changes of stratification in a deep water body with respect to its effect on phytoplankton sedimentation, a periodic function PIDI, with maxima in spring and fall and a minimum in mid-summer, is used. The geographical location of the water body is expressed not only by the water temperature, but also by the annual changes of the incident radiation (I) and photoperiodicity (Foxoa), giving the daily distribution of radiation. The following parameters characterize the individual water body: Q/V, ZNIX, TABLE 1 Model equations in the mathematical notation

a 1 = TEMP = 12+ 10 sin(t +220) a 2 = FOTOP =

12- cos(t)

a 3 = I = 280+ 210 sin(t + 240) a 4 = PIDI = 0.8 +

0.25 c o s ( t ) - 0.12 cos(2 t )

a 8 = 1.2a2/(blb2)

( 1.333 a r c t a n ( a 3 / ( 2 b 3 a 2 ) )

-- ( b 3 a 2 / a 3 ) ln(1

+

( a 3 / ( 2 b 3 a 2 ) ) 2)

- 1.333 arctan( a 3 exp( - bib 2 ) / ( 2 5 3 a 2 )) - b 3 a 2 / ( a 3 e x p ( - bib2) ) In(1 + ( a 3 e x p ( - b a b 2 ) / ( 2 6 3 a 2 ) ) 2) )0.366 exp(0.09al)

Phytoplankton £q = a s x l x 3 / (

Q o k 7 + x 3 ) - b s a l x 1 - b6XlX 2 - b4a4x I + v ( Y l - Xl)

Zooplankton Q

fC2 = b6/(lOOOk2)( k 4 k 8 X l X 2 ) / ( k 8 +

xl)_

x2(k

5 + k6)+

~ ( yO _

x2 )

Phosphorus Yc3 = k 3 a 8 x , x 3 / ( k +lO00k2x2(k

7 + x3)+ k3b6xlxz(1-

kn)ks/(k8

+ Xl)

Q o 5 + k 6 ) + v ( Y 3 - x3)

The meaning of the terms is as follows: a 1 temperature ( ° C), a 2 photoperiodicity (h), a 3 fight intensity (cal c m - 2 d a y - 1 ) , a 4 :sinking function (dimensionless), a 8 specific rate of integral phytoplankton photosynthesis (units of phosphorus per unit surface area per day). cal, calorie (International Table) = 4.1868 J (def).

9"7 TABLE 2 Parameter values for the standard run EPS ZMIX IK SINK

ba b2 b3 64

0.8 5

EFFUZ RZ

k4 k5

0.06 0.01

0.11

MORTZ

0.005

KSP

RESP

b5

KSF

0.075 200 1

FRZ

PRITA FRITZ

I,'10

Q~ V

b6 kI

0.005 0.01

k6 kv k8

0.05 0.01

K2 K3

kz k3

FRFOS

yO

0.03

0.05 1

yO

10

EFS extinction of water for light (units per m), EFFUZmaximum efficiency of phytoplankton utilization by zooplankton, ZMIXmixing depth (m), RZ specific rate of zooplankton respiration (per unit zooplankton and day), IK photosynthetic parameter (units of I), MORTZspecific zooplankton mortality rate (day-1), KSI"half-saturation constant for zooplankton grazing on algae (mg chlorophyU-a m-3), RESP specific phytoplankton respiration rate (day l)~ KSF half-saturation constant for phytoplankton uptake of phosphorus (mg P m-3), FRZ zooplankton grazing rate (m3 mg- 1 zooplankton day - 1), PRITAconcentration of phytoplankton in the inflow to the system, Q / V system throughflow rate (day a), FRITZ concentration of zooplankton in the inflow; FRrOS concentration of phosphorus in the inflow, K2 and K3 recalculation constants. EPS as well as the i n f l o w c o n c e n t r a t i o n s of p h o s p h o r u s (PRFOS), algae (PRITA) a n d z o o p l a n k t o n (PRITZ). T h e e q u a t i o n s of the m o d e l (1) are given in T a b l e 1. T a b l e 2 lists the p a r a m e t e r values used for the s t a n d a r d run. T h e p h o s p h o r u s flows into the w a t e r b o d y at a relative rate PRFOS ( Q / V ) a n d out of at a rate P ( Q / V ) . It is utilized b y p h o t o p l a n k t o n at a speed d e t e r m i n e d b y p h y t o p l a n k t o n g r o w t h rate a n d is released d u r i n g b o t h p h y t o - a n d z o o p l a n k t o n excretion. T h e p h y t o p l a n k t o n is c o n t r o l l e d b y light, t e m p e r a t u r e a n d p h o s p h o r u s , and its specific g r o w t h rate also d e p e n d s on p h y t o p l a n k t o n biomass. T h e b i o m a s s is grazed u p o n b y z o o p l a n k t o n a n d it sediments out of the w a t e r layer b e i n g c o n s i d e r e d ( p a r a m e t e r SINK). ZOOp l a n k t o n g r o w t h d e p e n d s o n the a m o u n t of p h y t o p l a n k t o n available (described b y the M i c h a e l i s - M e n t e n e q u a t i o n with p a r a m e t e r s FRZ and KSP). Z o o p l a n k t o n losses are d u e to m o r t a l i t y (MORTZ), b o t h f r o m n a t u r a l internal causes and fish grazing. GLOBAL BEHAVIOR FOR ZERO THROUGHFLOW

Constant conditions First we will e x a m i n e o u r s y s t e m G I R L u n d e r c o n s t a n t c o n d i t i o n s , i.e.

when TEMP, I, FOTOP a n d PIDI are c o n s t a n t .

98

Let us denote: V(x)

= x l q- x 2 -Jr-x 3 ,

V(x)

~ O

for

x>~O

Then follows: V ( x ( t ) ) = (1 - lO00k2)x2 [ ( b 6 k 4 / ( l O O O k 2 ) ) ( k S x l / ( k 8 + x l ) ) - k 5 - k6] -

xa(bsal + b4a,)

A natural condition, in the sense of the principle of conservation of matter, is that: 12(x) ~< 0

for all

x ~ R3+

To satisfy this condition it is sufficient to set: 1000k 2 = 1

(2)

The function 12(x) is a L y a p u n o v function on R3+ for the system (1). By the invariant principle we get that the limit set of each solution must lie in the maximal invariant subset M of the set (x ~ R 3 : l ) ( x ) = 0). In our case, M = (x ~ R3+: x~ = 0). When x ( t , x °) is a solution of the system (1), for t ~ ~ , x a ( t ) ~ O, system (1) on the set M has the following form: :~1 = 0 )c2 = - (ks + k 6 ) x 2 )c3 = (ks + k 6 ) x 2 ~o(x °) is invariant and therefore x 2 ~ 0 for t---, ~ , and in the aquatic system only phosphorus will be present. Let us denote S --- ( x ~ R 3, xa = x 2 = 0). The set S forms the set of all steady states of the system (1) for condition (2). Let us calculate the Jacobian of (1) in the steady state = (0, 0, c): a s c / ( k 7 + c ) - b s a 1 - b4a 4 A=

0

-a8c/(k7 + c)

- k5

-

k6

k5+ k6

Eigenvalues of the matrix A are as follows: v 1 = asc/(k

7 + c) - b,a 1 - b4a 4

v2 = _ k 5 - k 6

U3 ~-~-0

If v 1 > 0, i.e. a s c / ( k 7 + c) > b s a 1 + b4a 4, then the steady state is unstable. Biologically, this condition means that gains of the p h y t o p l a n k t o n are greater than its losses. In the opposite case the matrix A has two negative

99 eigenvalues and one zero value, then the multiple of zero eigenvalue is equal to the dimension of S and therefore there exists a neighbourhood U of 2, with the following characteristics: when x ° ~ U, then x ( t , x °) converges to the steady state from S. Seasonally' varying conditions

The functions TEMP, I, FOTOP and PIDI are periodical with a period 360, which approximates the effect of natural environment. From the definition of these functions it follows that they are bounded by positive constants. Under condition (2),

V(X(I)) <~ -Xl(b5almin q- b4a4min ) ~< 0 is valid. According to Yoshizawa (1975), x~(t) converges to zero for t ~ ~ . Therefore, for t sufficiently large, it holds that: k 2 = x2[ ( k 4 k 8 / ( k 8 + x , ) ) x , -

k 5 - k6] ~ - r x 2 <~ 0

By variation of constants we get x2(t ) --+ 0. GLOBAL BEHAVIOR FOR NONZERO THROUGHFLOW Constant conditions

Assume that Q / V > O and y O = 0 , i = 1 , 2, 3. In this case the point = (0, 0, 0) is the only steady state of (1). Let us denote z ( t ) = Xl(t) + Xz(t) + x3(t). By variation of constants we get:

where b = x ° + x ° + x ° and q = y O +yO2 +yO. For y,° = 0, z ( t ) - - - , O for all x ° ~ R 3, then x ( t , x °) ~ Y = (0, O, 0). The characteristic polynom of system (1) in the point ~ has the following negative eigenvalues: ~'1 = - b s a l

-

v2 = -k 5 -

k 6 --

Q/V

b4a4-

p/v

~:3 = - Q ~ v

The trivial steady state is assymptotically stable in the large. In the following we will demonstrate that for yO sufficiently small the existence of asymptotically stable, in the large, steady state remains. Let us denote by:

V(y)=(OFi(Y)) 9 Oxj

i.j=l

100 the Jacobian of the right hand side of (1) at y. There exists a 81 > 0, such that for all y ~ R 3, ]I Y I] <81, the matrix V ( y ) has eigenvalues with negative real parts. Let us consider the following inner product:

P(x, y)= (V(y) (x-y), x-y 5 It can be shown that there exists a 6 < 81 and a constant C > O, such that: P ( x , y) <~ (C - Q / V ) I I x - y II 2 for all x, y ~ S c l , where q v ~ - < 6 and S c = ( x : x l + x 2 + x 3~ i 0 ) The matrix V ( y ) is regular for II y II < 81, therefore Sc1 can contain a finite number of equilibria xa, . . . , x n only. Let us denote: Uxi=(x~U:x(t,x°)~xi

for

t~oo)

and

D=U-(U'x

1U, ..., UU'x,),

where U=intuS

c,

c3v~ < 8,

c 1 < c 3,

0
and U ;'x i C U x i

UXl is open, therefore there is a x ° ~ 8Ux I with the following properties: (1) X1 ~ O)(X 0) (2) x(t, x °) lies in D for t > 0. The set D does not contain equilibria. According to Pliss (1964), x ( t , x °) is stable and if z ° lies near to x °, then x(t, z °) converges to x ( t , x °) for t ~ oo. We can choose a z ° which is contained in Ux I as well. Therefore x a has to be asymptotically stable in the large. If S d does not contain equilibrium then there exist an asymptotically stable in the large periodic solution lying in Scl. The above argument is valid for all (c, Q / V ) , where 0 < c < c a and Q / V > C, evidently. We may now state the following: Proposition. Let Q~ V > C and CzV~-= 6. Suppose that 0 < c < c 2 and y is such that y~ > 0, i = 1, 2, 3, and c = yO + y20 + TO. Then the system (1) has an asymptotically stable in the large equilibrium or a nontrivial asymptotically stable in the large periodic solution. System (1) can be written in the form:

Yc= r ( x , yO) We have F(O, O) = 0 and the matrix:

- Q / V - bsa a - b4a 4 A =

0 0

0 - k 5 -- k 6 k 5 4- k

0 Q//V 6

0 -

Q/V

IUI

is regular. By the implicit function theorem there is a function g with the following properties: (1) g: U1 --+ U2, where U1 and U2 are neighbourhood of y0 = 0 and x = 0, respectively. (2) g(0) = 0, F(g(y°), y0) = 0 for all y0 ~ U1 Now we will demonstrate that the steady state Y = g ( y ° ) is positive for y,° small. Taylor expansion of the function g(y0) in point y0 = 0 has the following form: g ( y ° ) = By° + o( [I y° [I)

for

II y° LI --, 0

where Q/V Q / v + bsa 1 + b4a 4 B =

o

o

Q/v

0

k 5 ÷ k 6 4-

Q/V

(ks+k6)

0

Q/V(k s + k 6+ Q/V)

0 1

Evidently the following holds true: gi(Y°)>0

for

yO>0

and

yOsmall

The numerical results have shown that the steady state ~ exists for y o sufficiently large and yO > 0, and the Jacobian of (1) in point 2 has negative real parts, and therefore the existence of the asymptotically stable steady state remains. Variable conditions For periodic coefficients TEMP, I, FOTOP and PIDI under condition (4) we obtain that the solutions of (1) are ultimately bounded, i.e. there exist T in dependence of x ° that for all t > T, L[x(t, t o , x °) ]l < 3 q , where c I = v l ~ + ... +yO. According to Yoshizawa (1975) there exists a periodic solution with the period = 360. For details of the derivation see Kme] (1986). DEPENDENCE OF STEADY STATE ON PARAMETERS

The search of steady states was realized numerically by means of the computer EC 1045. Our effort was concentrated on the examination of the existence of stable steady states in dependence on the following parameters: (1) external

102

/

x

I t / t / I / iI r

• 5'0 . . . .

160. . . .

1~ . . . . 260. . . . 2~0. . . . 3~0. . . . 3~0. . . . *60. . . . . TIME

Fig. 1. Reaching steady state of equations (1) for external conditions corresponding to t = 210, starting from initial conditions Xl,0 = 100, x2,0 = 1, x3,0 = 20. Full line: phytoplankton concentration, x 1. Dashed line: zooplankton concentration, x 2. Dotted line: phosphorus concentration, x 3.

c o n d i t i o n s Q/V, ZMIX, EPS, (2) i n f l o w c o n c e n t r a t i o n o f p h o s p h o r u s , algae a n d z o o p l a n k t o n : PRFOS, PRITA a n d PRITZ, (3) i n t e r n a l s y s t e m p a r a m e t e r s : KSP, KSF, MORTZ. All were e x a m i n e d u n d e r c o n s t a n t c o n d i t i o n s c o r r e s p o n d ing to t = n × 30. A s t e a d y state for c o n s t a n t c o n d i t i o n s a n d a s e q u e n c e o f s t e a d y states are s h o w n in Figs. 1 a n d 2. T h e results for t = 210 which were selected as a p p r o x i m a t e l y r e p r e s e n t a t i v e o f the p e r i o d of m a x i m a l or m i n i m a l values of the variables xl, x2 a n d x 3 are d e m o n s t r a t e d in Figs. 3 to 5. F o r changes o f t h r o u g h f l o w conditions, Q/V (Fig. 3A) the relation is w i t h i n the r a n g e o f values b e t w e e n 0.001 a n d 0.1 p e r d a y ( c o r r e s p o n d i n g to r e t e n t i o n time o f w a t e r in the system 10 to 1000 days) r a p i d l y increasing for p h y t o p l a n k t o n , a n d the same is valid for p h o s p h o r u s . Z o o p l a n k t o n increases. System r e a c t i o n o n changes o f ZMIX is highly n o n l i n e a r (Fig. 3B). L o w effect is o b s e r v e d for 0 < ZMIX ~< 2, while a f t e r w a r d s p h y t o p l a n k t o n slightly increases a n d p h o s p h o r u s r a p i d l y decreases. Light c o n d i t i o n s as a f f e c t e d b y EPS show identical results (Fig. 3 C) T h e d e p e n d e n c e o n the p h o s p h o r u s i n f l o w c o n c e n t r a t i o n P R F O S studied within the r a n g e 1 0 - 2 0 0 m g p e r m 2 (Fig. 4A) shows a m o n o t o n o u s l y increasing u p w a r d c u r v e for p h y t o p l a n k t o n , w h e r e a s p h o s p h o r u s r a p i d l y

10 3

/ z

/ 6

. . . .

sb

. . . .

16o . . . .

1~o . . . .

~6o . . . .

(A)

2~o . . . .

36o . . . .

3~o

'

TIME

z

750

(B)

800

850

900

950

1000

1050

1100

TII~_

Fig. 2. Confrontation of a sequence of steady-state solutions for external conditions corresponding to successive time intervals and dynamic solution of equations (1). (A) Sequence of steady-state solutions. (B) Dynamic solution, periodic regime, See Fig. 1 for the meaning of the lines.

104

/

oloo

"

"

"

olo~

"

"

olo4

"

olo6

"

(A)

"

o:o8

o~

'

O/V

\ x

\ %

6

(B)

.

.

.

.

. 2.

.

~,

.

.

.

.

6.

.

.

.

.

8

tb

ZlIX

Fig. 3. Effect of external variables on the steady-state solution for external conditions corresponding to t = 210. (A) Effect of throughflow (Q/V) through the system. (B) Effect of the depth of the layer inhabited by algae (ZMIX). (C) Effect of the extinction coefficient of water for light, EPS. See Fig. 1 for the meaning of the lines. Due to scale identical for all variables, the changes of zooplankton are not seen.

105

\ \ \

\

j

f

J 0~0

(C)

o12

"

'

"

0'.4

'

"

"

o'.~

o.8

~.o

EPS

Fig. 3 (continued).

reaches a saturation level and zooplankton increases only slightly. The system does not seem to be sensitive on the values of inflow algal concentrations PRITA (investigated in the range 0.05 to 1 mg chlorophyll-a per litre) (Fig. 4B). On the other hand the inflow concentration of zooplankton (PRxxz, investigated in the range 0.001 to 0.1 mg C per litre) determines the steady state values much more: zooplankton steady-state concentration in the water body increases and surprisingly the same is valid for phytoplankton, while phosphorus barely changes (Fig. 4C). For the half-saturation constant of zooplankton feeding on phytoplankton (KsP), values within the range 0-200 mg 1-~ P were investigated. On day 210 the steady state values of phytoplankton react by a nonlinear decrease on increasing KSP, a near-saturation plateau being observed up to KSP of 50. Zooplankton is nearly constant, and phosphorus increases rapidly. The sensitivity on MORTZ is very low within the range 0-0.2 (Fig. 5C) DISCUSSION

We have studied the system GIRL under constant and naturally varying conditions without and with water inflow ( Q / V = 0 and Q/V > 0, respectively). As expected we have shown that, without water inflow, in the system GIRL only phosphorus remains and all living organisms vanish. The closed

106

40

30 ¸

20'

/

.

.

.

.

,.~

.

.

.

100

.

.

.

.

.

.

(A)

.

.

200

PRF'OS

-----

0.0

(B)

.150 .

o:2 .

.

.

. . . . . . . . . . . . . .

. 0.4. .

.

.

. o.e. .

.

.

. o.8.

~:o

PRITA

Fig. 4. Effect of the inflow concentrations on the steady state solution for external conditions corresponding to t = 210. (A) Effect of inflow phosphorus concentration, P~OS. (B) Effect of

the inflow concentration of phytoplankton, PRITA. (C) Effect of zooplankton concentration in the inflow, PRITZ. See Fig. 1 for the meaning of the lines.

107

J

~

J

J

- - ,

,

.

0.00

(c)

.

, 0.02

•

.

.

, 0.04

.

.

,

,

.

.

.

0.06

, 0.08

,

.

.

, 0.10

PRITZ

Fig. 4 (continued).

ecosystem is unable to exist, in spite of the energy supply to it. This is the same for both constant and variable coefficients. When water flows in (Q/V>O and y ° > 0 , i = 1, 2, 3), the system G I R L shows qualitatively different characteristics under constant and under naturally varying conditions. Under constant conditions the system reaches one steady state, which is for y/0 small asymptotically stable in the large. The numerical results have shown that this feature remains valid also for yO large. Under coefficient variation corresponding to natural conditions a steady state of this system does not exist. The effect of periodically varying environmental conditions produces a periodic trajectory with the 1-year periodicity. Figure 1 demonstrates the steady-state solution for constant conditions, and Fig. 2 shows a sequence of steady states calculated for every 30-day interval and a dynamic solution of the system. The comparison of Fig. 2A and B suggests that the trajectory of steadystate solutions is very close to the solution trajectory. However, they are not identical with the maximum of x 3. The behavior of the steady-state model solution (for t = 210) reactions on changes of input parameter values can be confronted with our knowledge about the behavior of natural aquatic ecosystems. The upward curve of phytoplankton biomass observed for the relation to PRFOS is in agreement with the Dillon and Rigler (1974) type of chlorophyll-a: total phosphorus

108

/ / / / / /

/

/

/

1

#

i

i i p

I

1 i

.

.

.

.

~

.

.

.

.

~

.

.

.

.

"~

(A)

KSP

(B)

KSF

.

.

.

.

250

Fig. 5. Effect of parameter changes on steady state solution of the system for conditions corresponding to t = 210. (A) Effect of the half-saturation constant for zooplankton feeding on phytoplankton, Ks]'. (B) Effect of the halfsaturation constant for phosphorus, KSF. (C) Effect of zooplankton mortality, MORTZ. See Fig. 1 for the meaning of the lines.

109

,

.

.

.

.

0.00

,

0.05

.

.

.

.

,

0.10

.

.

.

0.15

0.20

Fig. 5 (continued). relationships. However, to differ from both dynamic simulation results as well as natural observations evaluated in Stragkraba (1979), no saturation limit was observed. This is explainable while the concentrations of phosphorus in the water body at the retention time modelled (corresponding to about one year) are low. The behavior in respect to Q/V also is explainable in the light of limnological experiences. With increased phosphorus load connected with increasing Q relative to V, the phytoplankton biomass increases rapidly. We have already observed with dynamic simulations that although the increase of phosphorus load is due to both increased PRFOS and increased hydraulic load, the effect of both is not equivalent. The change of hydraulics also affects other system components directly, whereas phosphorus concentration acts only indirectly through increased phytoplankton growth. The phosphorus trajectory is smoothed out in the steady-state case and possesses a periodicity of nearly 2t. In total, the comparison suggests that the use of steady-state solutions for selected periods can help to study the system's reaction to various parameters. REFERENCES Dillon, P.J. and Rigler, F.H., 1974. The phosphorus-chlorophyll relationships in lakes. Limnol. Oceanogr., 19: 767-773.

110 Kme~, T., 1986. Dynamic ecological system models. Ph.D. thesis, Centre of Biological and Ecological Sciences, Slovak Academy of Sciences, Bratislava, 62 pp. (in Slovak). Pliss, V.A., 1964. Nonlocal Problems of Oscillation Theory. Nauka, Moscow, 205 pp. (in Russian). Stra~kraba, M., 1979. Mathematische Simulation der Produktionsdynamik in Gewassern und deren Anwendung auf die Produktionssteuerung in Talsperren. Z. Wasser Abwasser Forsch., 12: 56-64. Stra~kraba, M. and Gnauck, A., 1985. Freshwater Ecosystems. Modelling and Simulation. Developments in Environmental Modelling, 8. Elsevier, Amsterdam, 309 pp. Yoshizawa, T., 1975. Stability Theory and Periodic Solutions and Almost Periodic Solutions. Springer, New York, 233 pp.