Dynamics of the 'echo'effect in a phytoplankton system with nitrogen
Bulletin of Mathematical Biology 67 (2005) 487–507 www.elsevier.com/locate/ybulm
Dynamics of the ‘echo’ effect in a phytoplankton system with nitrogen fixation Khalid Boushaba∗,1, Mercedes Pascual Department of Ecology and Evolutionary Biology, University of Michigan, Ann Arbor, MI 48109-1048, USA Received 11 July 2003; accepted 17 August 2004
Abstract Consideration of nitrogen fixation adds a positive nonlinear feedback to plankton ecosystem models. We investigate the consequences of this feedback for secondary phytoplankton blooms and the response of phytoplankton dynamics to physical forcing. The dynamics of phytoplankton, Trichodesmium (the nitrogen fixer), and nutrients is modeled with a system of three differential equations. The model includes two types of nonlinear interactions: the competition of phytoplankton and Trichodesmium for light, and the positive feedback resulting from Trichodesmium recycling. A typical simulation of the model in time, with forcing by a varying mixed-layer depth, reveals a clear successional sequence including a secondary or ‘echo’ bloom of the phytoplankton. We explain this sequence of events through the stability analysis of three different steady states of the model. Our analysis shows the existence of a critical biological parameter, the ratio of normalized growth rates, determining the occurrence of ‘echo’ blooms and the specific sequence of events following a physical perturbation. The interplay of positive and negative feedbacks appears essential to the timing and the type of events following such a perturbation. © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved.
∗ Corresponding author.
E-mail address: [email protected] (K. Boushaba). 1 Present address: Department of Mathematics, 400 Carver Hall, Iowa State University, Ames, IA 50011, USA. 0092-8240/$30 © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.bulm.2004.08.004
488
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
1. Introduction Recent geochemical studies show deficiencies in N inputs relative to outputs in several ocean basins (Karl et al., 1992; Sambrotto et al., 1993; Bates et al., 1996, 1998; Michaels et al., 1996; Gruber and Sarmiento, 1997; Bates, 2001; Hood et al., 2001). New production estimates in the upper water column of tropical regions often exceed known new N inputs. Trichodesmium, the most prominent planktonic marine nitrogen fixer, is a major source of new nitrogen for the ocean (Capone and Carpenter, 1982; Capone et al., 1997, 1998; Capone, 2000; Karl et al., 2002). Nitrogen fixation is important ecologically, because nitrate is typically thought of as the limiting nutrient, but Trichodesmium use nitrogen rather than nitrate and nitrogen is in ample supply. Research through the 1970s indicated that Trichodesmium can be of regional significance but its contribution to the larger marine N and C cycle was considered minimal. Recent measurements from several different locations, including the tropical North Atlantic, the Arabian Sea and the Eastern Indian Ocean, have shown that N2 fixation by Trichodesmium can, at times, match or exceed the NO− 3 flux into the upper water column (MANTRA/PIRANA Cruise Data). Recently Hood et al. (2001) introduced an ecosystem model with six compartments to describe the effect of N2 on nitrogen and carbon fluxes at Bermuda Atlantic Time-Series Study (BATS). These compartments are no Trichodesmium phytoplankton, dissolved inorganic nitrogen, dissolved organic nitrogen, heterotrophs (zooplankton, microzooplankton, and bacteria), detritus, and Trichodesmium. In this model, Trichodesmium growth is controlled by the availability of light and its intrinsically slow growth rate. Trichodesmium interacts with no Trichodesmium phytoplankton by competing for light, which depresses its own growth, and by pumping nitrogen into the system, which stimulates no Trichodesmium phytoplankton growth. The ecosystem model has been embedded in a three dimensional circulation model of the Atlantic ocean from 45◦ N to 20◦ S (the Miami Isopycnal Coordinate Ocean or MICOM) (Coles et al., 2004a,b; Hood et al., in press). The coupled physical–biological model reproduces the general characteristics of the observed Trichodesmium biomass distribution. Furthermore, the model generates a clear successional sequence of events including a secondary phytoplankton bloom known as an ‘echo’ bloom. A recent analysis of remote sensing data provides support for the occurrence of such echo blooms fueled by nitrogen fixation (Coles et al., 2004a,b). A recent comparison of the coupled model with and without Trichodesmium suggests that it may be essential to explicitly represent both functional groups, the phytoplankton and the nitrogen fixers, in the ecosystem equations to accurately simulate the spatial and temporal variability of primary productivity and of new nitrogen inputs from nitrogen fixation (Coles et al., 2004a,b). To better understand the biological response of the system to physical forcing, we present in this paper an ecosystem model of reduced dimensionality capturing the essential components of the larger system by Hood et al. (2001). The model consists of three differential equations for the dynamics of phytoplankton, nitrogen fixers, and nutrients. It includes two types of nonlinear interactions between primary producers: the competition of phytoplankton and Trichodesmium for light, and the positive feedback resulting from Trichodesmium recycling.
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
489
Fig. 1. Model structure.
A typical simulation of the model in time, with forcing by a varying mixed-layer depth, reveals a clear successional sequence including a secondary no Trichodesmium phytoplankton bloom. We explain this sequence of events through the stability analysis of the three different steady states of the model. Our analysis shows the existence of a critical biological parameter, the ratio of normalized growth rates, which is non-dimensional and determinines the occurrence of ‘echo’ blooms and the specific sequence of events following a physical perturbation. 2. Model formulation The model consists of three coupled interacting components, describing changes in the concentrations of Nutrient (N) (mmol m−3 ), Phytoplankton (P) (mmol m−3 ) and Trichodesmium (T) (mmol m−3 ). The structure of the model is illustrated in Fig. 1, with arrows indicating flows of matter between the three compartments. The input to the system 0 is D Z (N − N) which represents nutrients entering from below the mixed layer. Arrows not ending in any compartment represent losses from the system. We assume that the mixed layer is thoroughly mixed at all times. Thus, the dynamics of the system can be represented by the three coupled differential equations, dP = U P P − SP P dt dT = G T T − ST T dt D dN = (N 0 − N) − U P P + ηS P P + γ ST T. dt Z We now discuss in detail the different terms and parameters of the model.
(1)
490
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
2.1. Phytoplankton growth The growth rate U P of phytoplankton is a function of light, nutrients and mixed layer depth. If light and nutrient limitation are independent, then we can write U P = I (z, P, T ) × µ P × u(N)
(2)
where I (z, P, T ) is the local light level, µ P is the maximum phytoplankton growth, I is the functional response of growth to light, and u(N) describes nutrient limitation. The local light level varies with z, the depth of the water column and the self-shading by phytoplankton (P), Thrichodesmium (T), and the attenuation by water itself. We then have (Steele, 1962) I = I (z, P, T ) = I0 e−[kx +k P ( P+T )]z
(3)
where I0 is the irradiance intensity at the sea surface, k P is the phytoplankton-specific attenuation coefficient, and k x is the short-wave extinction coefficient. The local light averaged over a mixed layer of depth Z , is given by Z ¯I = 1 I0 e−[kx +k P ( P+T )]z dz (4) Z 0 and I¯ =
I0 1 − e−[kx +k P ( P+T )]Z . Z [k x + k P (P + T )]
(5)
The function u(N) describes the nutrient uptake rate of phytoplankton and satisfies the following general assumptions: it is non-negative and increasing, saturates at high nutrient levels, and vanishes at zero (Hale and Somolinos, 1983). That is, u(N) is a continuous function defined on [0, ∞) satisfying u(0) = 0,
du >0 dN
and
lim u(N) = 1.
N→∞
(6)
In particular, this kind of function includes a Holling’s type II (Holling, 1959) or Michaelis–Menten function (Wroblewski and Richman, 1987) u(N) =
N , N+K
(7)
where K is the half-saturation constant. By neglecting the term e−[kx +k P ( P+T )]Z we write the production of new phytoplankton as N I0 µ P P × . Z [k x + k P (P + T )] N+K
(8)
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
491
2.2. Trichodesmium growth The growth of T is given by I0 µ T T (9) Z [k x + k P (P + T )] which includes the competition with P and other T cells for light through shading. It does not incorporate, however, any nutrient limitation because Trichodesmium fix nitrogen rather than utilize nitrate, and the nitrogen is plentiful in supply. 2.3. Loss and recycling of phytoplankton and Trichodesmium The loss of phytoplankton (Trichodesmium) is modeled by the expression S P P (ST T resp.), where S P and ST are the specific loss rates. These terms represent the mortality of cells due to death and grazing and the loss due to respiration. We assume that only a fraction of the dead phytoplankton (ηS P ) and the dead Trichodesmium (γ ST ) are recycled into dissolved nutrients. The model further assumes that the recycling is instantaneous, thus neglecting the time required to regenerate nutrient from dead biomass by bacterial decomposition. However, a delay in nutrient recycling is present in natural systems and varies inversely with temperature (Whitttaker, 1975). Model (1) can easily be modified to incorporate these delays in the nutrient equation, t D dN = (N 0 − N) − U P P + ηS P f (t − s)P(s)ds dt Z −∞ t g(t − s)T (s)ds. (10) + γ ST −∞
The delay-kernels f (s) and g(s) are non-negative bounded functions defined on [0, ∞) and describe the contribution to the nutrient recycled at time t of phytoplankton and Trichodesmium cells that died in the past. We normalize the kernels such that ∞ ∞ f (s)ds = 1, g(s)ds = 1. (11) 0
0
We note that the presence of the distributed time delays does not affect the equilibrium values. For simplicity, we will assume that the nutrient recycling of the phytoplankton and Trichodesmium is instantaneous, so that f (s) = g(s) = δ(s), the Dirac delta function. 2.4. Mixed-layer depth We assume that the mixed-layer depth Z varies during the year as follows, 0 ≤ t < 60 200 60 ≤ t ≤ 270 Z = 20 2t − 520 270 ≤ t ≤ 360, where t is measured in days and Z , in meters. In the simulations, we consider that the mixed layer-depth varies continuously throughout the year, except for a discontinuous jump of 180 m at t = 60. This piecewise-linear function provides a simple approximation of the three major features of the annual cycle (Evans and Parslow, 1985; Fasham, 1995; Hood
492
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
et al., 2001): a rapid shallowing in the spring, a constant shallow layer throughout the summer, and a deepening during the autumn and winter. 2.5. Mixing The model assumes that the ocean is divided into two completely mixed layers: an upper one, biologically active, which contains nutrients and plankton of constant concentration 0 with depth, and a lower one, inactive, with nutrients but no plankton. The D Z (N − N) term models the exchange of nutrients with the water below the mixed layer, which is assumed to have a constant concentration N 0 . The mixing occurs at a constant rate D (m day−1 ). D Z represents the fraction of the mixed layer that is exchanged with deep water through diffusive processes. 3. Interpretation of the ‘echo effect’ Before proceeding to the mathematical analysis of the model, we illustrate with three representative simulations the biological response of the system to changes of the mixedlayer depth. These simulations correspond to different values of the non-dimensional critical parameter ζ =:
µT ST µP SP
relating the maximum growth rate of T to that of P. The growth rates in ζ are normalized by the respective loss rates of these two groups. We show here that the different responses in the simulations can be explained through the existence of three different steady states whose local stability properties vary with ζ . The existence and stability properties of these steady states are studied in detail in the next section (Section 4). Fig. 2(a) shows a typical simulation of the system for a low growth rate of T relative 0 to that of P (ζ < N 0N+K < 1). For this choice of parameters, we observe the typical successional sequence already described for the higher dimensional model (See Table 1). In this sequence, mixed-layer depth shallowing stimulates a phytoplankton bloom, which is in turn followed by a Trichodesmium bloom after nutrients become depleted. The lack of nutrients and the concomitant decrease of P, releases T from competition allowing it to grow in spite of its maximum low growth rate. The Trichodesmium bloom then stimulates phytoplankton growth and elevates phytoplankton biomass. We refer here to this specific successional pattern as the ‘echo effect’. Fig. 2(b) shows the corresponding phase-space diagram and the three steady states E 0 = (0, 0, N 0 ), E 2 = (P2 , 0, N2 ), and E 3 = (P ∗ , T ∗ , N ∗ ). As the successional pattern proceeds, the trajectory of the system moves from the proximity of E 0 , to that of E 2 , and finally, E 3 . In the winter, the values of P and T remain low and the nutrient concentration approaches that of the sub-mixed-layer water, N 0 . We will show in Section 4 that the steady state E 0 = (0, 0, N 0 ) is globally stable provided the mixed layer is sufficiently deep, as is indeed the case in the winter. In the spring, a bloom of P occurs, as T remains low and N is almost entirely consumed. We will show that, as the mixed layer shallows, the steady state E 2 = (P2 , 0, N2 ) becomes locally stable. For the even shallower summer mixed layer,
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
493
Fig. 2. (a) Temporal dynamics of the model showing an ‘echo’ bloom: phytoplankton (P) (dashed line), Trichodesmium (T) (dashed–dotted line), nutrients (N) (thick line). The mixed-layer depth MLD is in meters (solid line). (b) Phase-space diagram.
E 2 loses stability and E 3 = (P ∗ , T ∗ , N ∗ ) becomes in turn locally stable. The tendency of the system towards E 3 explains the bloom of T and the secondary bloom of P. Finally, by the end of the summer, the deepening of the mixed layer returns the system towards E 0 . Fig. 3 shows the simulated time series for a larger growth rate of Trichodesmium such that N0 < ζ < 1, N0 + K with all other parameters kept at their previous values. In this case, the trajectory in phase space goes from E 0 directly to E 1 = (0, T1 , N1 ), then to E 3 , returning to E 0 . The absence of E 2 which remains unstable at all times, underlies the lack of a spring bloom of P. The system responds to the shallowing of the mixed layer directly with a bloom of T and the echo effect does not occur. There is instead a single P bloom stimulated in part by nitrogen fixation. Fig. 4 shows the simulated times series for an even larger growth rate of Trichodesmium, such that ζ > 1. In this case Trichodesmium blooms in the spring but both phytoplankton blooms are absent. Biologically the phytoplankton spring bloom does not occur because
494
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
Fig. 3. (a) Temporal dynamics with a single phytoplankton bloom following that of the nitrogen fixers. There is no phytoplankton spring bloom. (b) Phase-space diagram (see text).
the Trichodesmium bloom first and then the phytoplankton become light limited, as the result of the negative feedback between these groups. The trajectory in phase space goes from the proximity of E 0 to E 1 = (0, T1 , N1 ), returning to E 0 . The above cases, as well as a more general picture of the system’s dynamics, are summarized in a two-dimensional bifurcation diagram (Fig. 5). This diagram illustrates the stability of the equilibria as a function of two parameters: mixed layer depth, MLD, and the ratio of maximum growth rates ζ . The functions and parameters used to plot the different regions in this graph are given in the next section. To predict the specific response to a shallowing of the mixed layer, a vertical line can be drawn for the relevant value of ζ . The regions crossed by this line tell us whether an echo effect will occur, or alternatively, whether the main bloom of P will be absent. As we will show, however, this graph gives us the global dynamics of E 0 , E 1 , and E 2 but not that of the interior steady state E 3 . To complete the picture, we need to consider the additional parameter µT . Fig. 6 shows that for ζ = 0.1374, the mixed layer fixed at 20 m and all other parameters as in the simulation of Fig. 2 except for an extremely large (and biologically unrealistic) value of µT = 45.9, the trajectory in phase space goes to a periodic orbit shown in Fig. 6. This result, which is relevant to the appearance of multiple blooms, will be proven mathematically at the end of next section.
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
495
Fig. 4. (a) Temporal dynamics with no phytoplankton bloom following the shallowing of the mixed-layer depth. (b) Phase-space diagram.
4. Steady-state analysis System (1) takes the specific form I0 µ P P N dP = − SP P dt Z [k x + k P (P + T )] N + K dT I0 µ T T = − ST T dt Z [k x + k P (P + T )] dN D N I0 µ P P = (N 0 − N) − + ηS P P + γ ST T, dt Z Z [k x + k P (P + T )] N + K
(12)
with initial conditions P(0) = P0 ≥ 0,
T (0) = T0 ≥ 0,
N(0) = N0 ≥ 0.
By adding all the three equations, we can see that d (P + T + N) = −D N − (1 − η)S P P − (1 − γ )ST T dt D I0 µ T T + N0 + Z [k x + k P (P + T )] Z
(13)
496
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
Fig. 5. Bifurcation diagram as a function of mixed-layer depth (MLD) and normalized growth rate ζ (see text).
≤ −D1 (P + T + N) +
I0 µ T D + N0 ZkP Z
where D1 = max D Z , (1 − η)S P , (1 − γ )ST . Thus, the solutions are bounded. From the first equation of system (12), it also follows that P = 0 is an invariant subset, that is, P = 0 if and only if P(t) = 0 for some t. Thus, P(t) > 0. A similar argument holds for T = 0 using the second equation of system (12). We further have N(t) positive. Indeed, if N(t) = 0 for some t, then dN dt > 0 which implies that N remains positive. This proves: Proposition 1. All solutions of system (12) under (13) are bounded and positive. We begin with equilibrium E 0 = (0, 0, N 0 ) which exists for all parameter values. The local stability of E 0 is determined by the eigenvalues of the Jacobian matrix evaluated at this point, I0 µ P N 0 0 0 − SP Z kx N 0 + K I µ 0 T M0 = 0 − ST 0 . Z kx I0 µ P N 0 − γ S −D + ηS P T 0 Z kx N + K The eigenvalues of M0 are I0 µ P N 0 I0 µ T D y2 = − ST , y3 = − . − SP , 0 Z kx N + K Z kx Z Hence, the local stability of E 0 depends on the depth of the mixed layer as described below. This proves: y1 =
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
497
Fig. 6. (a) Temporal dynamics showing periodic solutions. The concentration of T is magnified (×20). (b) Phasespace diagram.
Proposition 2. If I0 µ T Z > max ψ(ζ ), , ST k x
(14)
where ψ(ζ ) =
I0 µ T N0 1 ST k x N 0 + K ζ
then E 0 is asymptotically stable. By changing slightly condition (14) to I0 µ P I0 µ T Z > max , , S P k x ST k x
(15)
we can establish the global stability of E 0 using an a priori estimate of P, T , and N. First we can see that dP I0 µ P P N = − SP P dt Z [k x + k P (P + T )] N + K
498
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
≤
I0 µ P − S P P. Z kx
which using Gronwall’s inequality (Hale, 1969, p. 35), gives I0 µ P P(t) ≤ P0 exp − SP t . Z kx I0 µ P S P kx ,
limt →∞ P(t) = 0. A similar argument follows for T with I0 µ T − ST t . T (t) ≤ T0 exp Z kx
Thus, if Z >
Thus, if Z >
I0 µT S T kx
(16)
, limt →∞ T (t) = 0. On one hand we have
I0 µ P P D N dN = (N 0 − N) − + ηS P P + γ ST T dt Z Z [k x + k P (P + T )] N + K D ≤ (N 0 − N) + ηS P P + γ ST T Z D 0 ≤ (N − N) + (ηS P P0 + γ ST T0 ) exp(−δt) Z
(17)
I0 µ T I0 µ P − SP , − − ST . δ = min − Z kx Z kx
(19)
(18)
where
On the other hand, we also have D N I0 µ P P dN = (N 0 − N) − + ηS P P + γ ST T dt Z Z [k x + k P (P + T )] N + K D I0 µ P ≥ (N 0 − N) − P. Z kx Z
(20)
Let h = N − N 0 . Then, by combining inequalities (18), (20) and using similar arguments to those already used for P, it follows that limt →∞ h(t) = 0. In conclusion, under condition (15), limt →∞ (P(t), T (t), N(t)) = E 0 . This proves: Proposition 3. Under condition (15) E 0 is globally stable. By setting P = 0, ddtT = 0 and for equilibrium E 1 ,
dN dt
I0 µ T = ST Z (k x + k P T ) D 0 (N − N) + γ ST T = 0. Z The solutions are given by T1 =
I0 µ T kx − Z k P ST kP
= 0, we next obtain the following system of equations
(21) (22)
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
N1 = N 0 +
γ ST
D Z
I0 µ T kx − Z k P ST kP
499
.
The steady state N1 is interesting because N1 > N 0 , unlike the case for normal NPZ models (Edwards, 1997; Edwards and Brindley, 1997). Biologically, this difference shows that the presence of Trichodesmium increases the nitrate steady-state value above that for which Trichodesmium is absent, because of nitrogen fixation. The equilibrium E 1 = (0, T1 , N1 ) is positive if and only if Z < kIx0SµTT . Notice that this condition immediately implies the instability of E 0 . Furthermore, as Z approaches kIx0SµTT , E 1 approaches E 0 . The Jacobian matrix at E 1 is µ P ST N1 − SP 0 0 µT N1 + K I0 µ T Z ST k x kPT M1 = − − 1 ST 0 2 Z (k x + k P (P + T )) I0 µ T I0 µ P kx + k P T N D + ηS P − γ ST − Z (k x + k P (P + T ))2 N + K Z with eigenvalues N1 µ P ST y1 = − SP , µT N1 + K
y2 =
Z ST k x − 1 ST , I0 µ T
By writing µT S P , µ P ST
ζ =
it follows that a necessary condition for the stability of E 1 is N0 . N0 + K This proves: ζ >
Proposition 4. If one of the following conditions is satisfied (i)
N0 N 0 +K
< ζ < 1 and ϕ(ζ ) < Z
1
and
Z
N0 , N 0 +K
E 2 is locally asymptotically stable if and only if I0 µ T D ζK − N0 Z< + and ζ > 1. k x ST λk x 1 − ζ
Hereafter we assume that ζ < 1. The interior steady state E 3 = (P ∗ , T ∗ , N ∗ ) is then given by Kζ 1−ζ γ ST I0 µT D 1 P∗ = − k x + (N 0 − N ∗ ) λ + γ ST k P Z ST Z 1 λ I D µ 0 T − k x − (N 0 − N ∗ ) . T∗ = λ + γ ST k P Z ST Z N∗ =
Furthermore the following condition Z < min(ϕ(ζ ), φ(ζ )) ensures the positivity of (P ∗ , T ∗ ), with I0 µ T ζK Dk P − N0 φ(ζ ) = + . ST k x 1−ζ λk x ST The Jacobian matrix at E 3 is given by m 1 m 1 m 13 M = m2 m2 0 m 31 m 32 m 33 where
Z m 1 = −k P S P ST P∗ < 0 I0 µ T Z m 2 = −k P ST2 T∗ < 0 I0 µ T µP K P∗ m 13 = ST >0 µT (K + N ∗ )2
(27) (28) (29)
502
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
m 31 = −(k x + k P T ∗ )ST S P
Z + ηS P I0 µ T
m 32 = −m 1 + γ ST > 0 D m 33 = − − m 13 < 0. Z The corresponding characteristic equation is λ3 + a1 λ2 + a2 λ + a3 = 0
(30)
where a1 = −(m 1 + m 2 + m 33 ) a2 = m 33 (m 1 + m 2 ) − m 31 m 13 a3 = m 2 m 13 (m 31 − m 32 ). By the Routh–Hurwitz criteria, all roots of Eq. (30) have negative real parts if and only if a1 > 0,
a3 > 0
and
a1 a2 − a3 > 0.
We have that a1 > 0 and a3 = m 2 (η − 1) > 0. We also have a1 a2 − a3 = −(m 1 + m 2 + m 33 )(m 33 (m 1 + m 2 ) − m 31 m 13 ) − m 2 (η − 1) = −(m 1 + m 2 + m 33 )m 33 (m 1 + m 2 ) + (m 1 + m 2 + m 33 )m 31 m 13 − m 2 (η − 1) = −(m 1 + m 2 + m 33 )m 33 (m 1 + m 2 ) + (m 1 + m 33 )m 31 m 13 + m 2 (m 31 m 13 − (η − 1)). Thus, if m 31 < 0
and
(m 31 m 13 − (η − 1)) < 0
then a1 a2 − a3 > 0. One can see that if Z > θ (ζ ) then (m 31 m 13 − (η − 1)) < 0, where θ (ζ ) =
I0 µ T k x ST
γ ST (1−η) 2 k P (1 − ζ ) − S P ST (λ + γ ST )K γ ST D 2 0 k P (1 − ζ ) + Z (N + K )
.
The function θ is negative for the model parameters chosen for the main run (Fig. 2). This proves: Proposition 6. The positive steady state E 3 is locally stable if η
I0 µ T < Z < min(ϕ(ζ ), φ(ζ )). k x ST
Now assume that the steady state E 3 is locally stable. We would like to know if E 3 loses its stability when one of the parameters is varied. Note the local stability of E 3 is equivalent to the conditions a1 > 0, a2 > 0 and a1 a2 > a3 . Thus, we would like to know whether
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
503
these conditions will be violated when one of the parameters changes. More specifically, we would like to determine if the real part of the roots of the characteristic equation (30) changes from negative to zero to positive, that is, if a Hopf bifurcation occurs at E 3 . We choose µT , the maximal growth rate of Trichodesmium, as the bifurcation parameter, and address the existence of a critical value µ¯ T such that a1 (µ¯ T )a2 (µ¯ T ) = a3 (µ¯ T ), d [a1 (µT )a2 (µT ) − a3 (µT )] < 0. dµT µT =µ¯ T The characteristic equation (30) becomes (λ + a1 (µ¯ T ))(λ2 + a2 (µ¯ T )) = 0, which has a negative real root λ1 = −a1 (µ¯ T ) and a pair of purely imaginary roots λ2,3 = ±i a2 (µ¯ T ). By continuity, for µT ∈ (µ¯ T − θ, µ¯ T + θ ), with θ > 0 sufficiently small, the eigenvalues are of the form λ1 = −a1 (µT );
λ2,3 = σ (µT ) ± iν(µT ).
To determine if a Hopf bifurcation occurs when µT = µ¯ T , we need to verify the transversality condition. By substituting λ2 into the characteristic equation (30), calculating the derivative with respect to µT , and separating the real and imaginary parts, we obtain dσ (µT ) dν(µT ) − B(µT ) + C(µT ) = 0, dµT dµT dσ (µT ) dν(µT ) B(µT ) + A(µT ) + R(µT ) = 0, dµT dµT A(µT )
where A(µT ) = 3(σ 2 (µT ) − ν 2 (µT )) + 2a1(µT )σ (µT ) + a2 (µT ), B(µT ) = 6σ (µT )ν(µT ) + 2a1 (µT )ν(µT ), C(µT ) = a1 (µT )(σ 2 (µT ) − ν 2 (µT )) + a22 (µT )σ (µT ) + a3 (µT ), R(µT ) = 2a1 (µT )σ (µT )ν(µT ) + a2 (µT )ν(µT ).
Notice that σ (µ¯ T ) = 0, ν(µ¯ T ) = have dRe λ2,3 = σ (µ¯T ) dµ T
µT =µ¯T
= − > 0.
d dµT
√ aµ¯ T . By solving for
dσ dµT
[a1 (µT )a2 (µT ) − a3 (µT )] 4[a12(µT ) + a2 (µT )]
and setting µT = µ¯ T , we
µT =µ¯T
504
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
Table 1 Model parameter symbols, units and values used in the simulations of the echo bloom Description Biological parameters Maximum phytoplankton growth rate Maximum Trichodesmium growth rate Phytoplankton natural mortality rate Trichodesmium natural mortality rate Sat. const. for N uptake by P Phytoplankton-specific attenuation coeff. Recycling efficiency of phytoplankton Recycling efficiency of Trichodesmium Light efficiency at the sea surface Physical parameters Shortwave extinction coefficient Concentration below mixed layer Mixing rate Mixed layer depth
Symbol
Value
Units
µP µT SP ST K kP η γ I0
3.22 0.17 0.08 0.09 0.5 0.0223 0.5 0.8 0.116
day−1 day−1 day−1 day−1 mmol m−3 m2 mmol−1 Dimensionless Dimensionless Dimensionless
kx N0 D Z
0.03 0.6 0.005 20–200
m−1 mmol m−3 day−1 m
Thus, the real part of λ1,2 changes from negative to zero at µT = µ¯ T and becomes positive when µT > µ¯ T . A summary of the above analysis is given in the following proposition on the bifurcation of system (12). Proposition 7. If there exists a positive number µ¯ T such that a1 (µ¯ T )a2 (µ¯ T ) = a3 (µ¯ T )
and
d [a1 (µT )a2 (µT ) − a3 (µT )] µ¯ T , E 3 becomes unstable and a family of periodic solutions bifurcates from E 3 . 5. Discussion We have presented a model with two functional groups, the traditional ‘compartment’ P and an additional one representing nitrogen fixers. The model includes two types of nonlinear interactions between these groups: the competition of phytoplankton and Trichodesmium for light, and the positive feedback resulting from nutrient inputs via Trichodesmium. Although quite simple, this system captures the essential features of the larger ecosystem model introduced by Hood et al. (2001) and used in coupled biological–physical simulations (Hood et al., in press). Through the stability analysis of three different equilibria, we have shown the existence of a critical biological parameter ζ , the ratio of normalized growth rates. A two-dimensional bifurcation diagram based on ζ and mixedlayer depth, was developed to explain the different possible paths of the system in response to the yearly seasonal variation of mixed-layer depth. Mixed-layer depth has long been recognized as a critical physical parameter for the occurrence of phytoplankton blooms. Our
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
505
analysis shows that, when nitrogen fixation is considered, the biological parameter ζ is also important, specifically to the occurrence of a secondary ‘echo’ bloom. As for other models, the first bloom in the spring occurs in response to the physics. The second one is then due to the biological interactions and occurs when ζ is small. It is interesting to note that this second bloom occurs in the autumn. Thus ‘echo’ blooms fueled by nitrogen fixation provide an alternative explanation for phytoplankton increases observed at that time of the year, typically explained by the deepening of the mixed layer while sunlight is still plentiful. A third axis not included in the bifurcation diagram concerns the maximum growth rate of the nitrogen fixers. For fast growth, it is still possible to retain a low ζ provided the mortality rate is also high. Under these conditions, we have shown the tendency of the system to oscillate. Even decaying oscillations would set a scenario in which multiple phytoplankton blooms follow the main spring bloom. Although parameters in this region appear so far unrealistic, the feasibility of this scenario remains open with the recent discovery of small nitrogen fixers in the picoplankton component, whose identity and biology is only beginning to be understood (Zehr et al., 2001). Two potentially important limiting nutrients are not represented in our model: phosphorus and iron. Fe is required for nitrogenase, the enzyme responsible for N2 reduction activity. Because Fe is extremely insoluble in oxygenated seawater, the bioavailability of Fe is controlled by the presence of a variety of organic ligands that enhance Fe solubility (Rue and Bruland, 1995; Karl et al., 2002). Questions therefore arise concerning the proximate control(s) of N2 -fixation rates, including Fe, phosphorus, and other growth factors. Although such controls have been omitted from our equations, our results and general approach may prove useful in interpreting the output of alternative, more complicated models. In particular, it will be interesting to ascertain whether a similar bifurcation diagram underlies the trajectories of alternative models in which a factor other than light controls growth and underlies the competitive interaction of the two functional groups. For environments in which competition for light is not critical (e.g. the Pacific ocean, T. Michaels, personal communication), we propose that the mixed-layer depth axis in the diagram would be replaced by the new limiting factor (such as iron). This possibility remains to be determined. Another extension concerns the treatment of explicit space and the relevant spatial scales at which to address the spatio-temporal dynamics of phytoplankton ecosystems with nitrogen fixation. The large numerical simulations of coupled biological–physical models (Coles et al., 2004a; Hood et al., in press) have not so far resolved space sufficiently to address the effect of eddies, the mesoscale variability or ‘weather’ of the ocean. This question is particularly intriguing given the recent discovery of small nitrogen fixers (Fuhrman and Capone, 2001; Zehr et al., 2001) with potentially faster growth rates, and therefore also faster temporal scales of response to physical forcing. Little is known about their role in the dynamics of plankton ecosystems. Acknowledgements We thank Raleigh Hood for early discussion on the model formulation. This work was partially supported by an NSF Biocomplexity grant.
506
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
References Bates, N.R., 2001. Interannual variability of oceanic CO2 and biogeochemical properties in the Western North Atlantic subtropical gyre. Deep-Sea Res. II 48, 1507–1528. Bates, N.R., Micheals, A.F., Knap, A.H., 1996. Seasonal and interannual variability of oceanic carbon dioxide species at the US JGOFS Bermuda Atlantic Time-series Study (BATS) site. Deep-Sea Res. 43, 347–383. Bates, N.R., Takahashi, T., Chipman, D.W., Knap, A.H., 1998. Variability of pCO2 on diel to seasonal timescales in the Sargasso Sea near Bermuda. J. Geophys. Res. 96, 15567–15585. Capone, D.G., 2000. The marine nitrogen cycle. In: Kirchman, D. (Ed.), Marine Microbial Ecology. John Wiley, NY, pp. 455–493. Capone, D.G., Carpenter, 1982. Nitrogen fixation in the marine environment. Science 217, 1140–1142. Capone, D.G., Subramaniam, A., Montoya, J.P., Voss, M., Humborg, C., Johansen, A.M., Siefert, R.L., Carpenter, E.J., 1998. An extensive bloom of the N2-fixing cyanobacterium, Trichodesmium erythraem in the central Arabian Sea. Mar. Ecol. Prog. Ser. 172, 281–292. Capone, D.G., Zehr, J.P., Pearl, H.W., Bergman, B., Carpenter, E.J., 1997. Trichodesmium, a globally significant marine cyanobacterium. Science 276, 1221–1229. Coles, V.J., Hood, R.R., Pascual, M., Capone, D.G., 2004a. Modeling the impact of Trichodesmium and nitrogen fixation in the Atlantic ocean. Ocean. J. Geophys. Res. 109, C06007. Coles, V.J., Wilson, C., Hood, R.R., 2004b. Remote sensing of new production fuelled by nitrogen fixation. Geophys. Res. Lett. 31, L06301. Edwards, A.M., 1997. A rational dynamical-systems approach to plankton population modeling. Ph.D. Thesis, University of Leeds. Edwards, A.M., Brindley, J., 1997. Oscillatory behavior in a three-component plankton population model. Dyn. Stability System 11, 347–370. Evans, G., Parslow, J., 1985. A model of annual plankton cycles. Biol. Oceanogr. 3, 327–347. Fasham, M.J.R., 1995. Variations in the seasonal cycle of biological production in subarctic oceans: a model sensitivity analysis. Deep-Sea Res. 42, 1111–1149. Fuhrman, J.A., Capone, D.G., 2001. Nifty nanoplankton. Nature 412, 593–594. Gruber, N., Sarmiento, J.L., 1997. Global pattern in marine nitrogen fixation and denitrification. Global Biogeochem. Cycles 11, 235–266. Hale, J.K., 1969. Ordinary Differential Equations. Wiley-interscience, New York. Hale, J.K., Somolinos, A.S., 1983. Competition for fluctuating nutrients. J. Math. Biol. 18, 255–280. Holling, C.S., 1959. Some characteristics of simple types of predation and parasitism. Canad. Entomol. 91, 385–398. Hood, R.R., Capone, D.G., Olson, D.B., 2001. Modeling seasonal and interanual biogeochemical and N2 cycles at BATS. Deep-Sea Res. II 48, 1609–1648. Hood, R.R., Coles, V.J., Capone, D.G., 2004. Modeling the Distribution of Trichodesmium and Nitrogen Fixation in the Atlantic, Ocean. J. Geophys. Res. (in press). Karl, D. et al., 2002. Dinitrogen fixation in the world’s oceans. Biogeochemistry 57, 517–519. Karl, D., Letelier, R., Hebel, D.V., Bird, D.F., Winn, C.D., 1992. Trichodesmium blooms and new nitrogen in the North Pacific gyre. In: Capone, D.G., Reuter, J.G. (Eds.), Marine Pelagic Cyanobacteria: Trichodesmium and other Diazotrophs. Kluwer Academic Publishers, Boston, pp. 219–237. Michaels, A.F., Olson, D., Sarmiento, J., Ammerman, J., Fanning, K., Jahnke, R., Knap, A.H., Lipschultz, F., Prospero, J., 1996. Inputs, losses and transformations of nitrogen and phosphorus in the Pelagic North Atlantic Ocean. Biogeochemistry 35, 181–226. Rue, E., Bruland, K., 1995. Complexation of iron (III) by natural organic ligands in the Central North Pacific as determineted by new competitive ligand equilibration/adsorptive cathodic voltametric method. Mar. Chem. 117–138. Sambrotto, R.N., Savidge, G., Robinson, C., Boyd, P., Takahshi, T., Karl, D.M., Langdon, C., Chipman, D., Marra, J., Codispoti, L., 1993. Elevated consumption of carbon relative to nitrogen in the surface ocean. Nature 363, 248–250. Steele, J.H., 1962. Environmental control of photosynthesis in the sea. Limnol. Oceanogr. 7, 137–150. Whitttaker, R.H., 1975. Communities and Ecosystems. Macmillan, New York.
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
507
Wroblewski, J.S., Richman, J.G., 1987. Nonlinear response of plankton to wind mixing events—implications for survival of larval northern anchovy. J. Plank. Res. 9, 103–123. Zehr, J.P., Waterbury, J.B., Turner, P.J., Montoya, J.P., Omoregie, E., Steward, G.F., Hansen, A., Karl, D.M., 2001. Unicellular cyanobacteria fix N2 in the subtropical North Pacific Ocean. Nature 412, 635–638.
Dynamics of the ‘echo’ effect in a phytoplankton system with nitrogen fixation Khalid Boushaba∗,1, Mercedes Pascual Department of Ecology and Evolutionary Biology, University of Michigan, Ann Arbor, MI 48109-1048, USA Received 11 July 2003; accepted 17 August 2004
Abstract Consideration of nitrogen fixation adds a positive nonlinear feedback to plankton ecosystem models. We investigate the consequences of this feedback for secondary phytoplankton blooms and the response of phytoplankton dynamics to physical forcing. The dynamics of phytoplankton, Trichodesmium (the nitrogen fixer), and nutrients is modeled with a system of three differential equations. The model includes two types of nonlinear interactions: the competition of phytoplankton and Trichodesmium for light, and the positive feedback resulting from Trichodesmium recycling. A typical simulation of the model in time, with forcing by a varying mixed-layer depth, reveals a clear successional sequence including a secondary or ‘echo’ bloom of the phytoplankton. We explain this sequence of events through the stability analysis of three different steady states of the model. Our analysis shows the existence of a critical biological parameter, the ratio of normalized growth rates, determining the occurrence of ‘echo’ blooms and the specific sequence of events following a physical perturbation. The interplay of positive and negative feedbacks appears essential to the timing and the type of events following such a perturbation. © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved.
∗ Corresponding author.
E-mail address: [email protected] (K. Boushaba). 1 Present address: Department of Mathematics, 400 Carver Hall, Iowa State University, Ames, IA 50011, USA. 0092-8240/$30 © 2004 Society for Mathematical Biology. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.bulm.2004.08.004
488
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
1. Introduction Recent geochemical studies show deficiencies in N inputs relative to outputs in several ocean basins (Karl et al., 1992; Sambrotto et al., 1993; Bates et al., 1996, 1998; Michaels et al., 1996; Gruber and Sarmiento, 1997; Bates, 2001; Hood et al., 2001). New production estimates in the upper water column of tropical regions often exceed known new N inputs. Trichodesmium, the most prominent planktonic marine nitrogen fixer, is a major source of new nitrogen for the ocean (Capone and Carpenter, 1982; Capone et al., 1997, 1998; Capone, 2000; Karl et al., 2002). Nitrogen fixation is important ecologically, because nitrate is typically thought of as the limiting nutrient, but Trichodesmium use nitrogen rather than nitrate and nitrogen is in ample supply. Research through the 1970s indicated that Trichodesmium can be of regional significance but its contribution to the larger marine N and C cycle was considered minimal. Recent measurements from several different locations, including the tropical North Atlantic, the Arabian Sea and the Eastern Indian Ocean, have shown that N2 fixation by Trichodesmium can, at times, match or exceed the NO− 3 flux into the upper water column (MANTRA/PIRANA Cruise Data). Recently Hood et al. (2001) introduced an ecosystem model with six compartments to describe the effect of N2 on nitrogen and carbon fluxes at Bermuda Atlantic Time-Series Study (BATS). These compartments are no Trichodesmium phytoplankton, dissolved inorganic nitrogen, dissolved organic nitrogen, heterotrophs (zooplankton, microzooplankton, and bacteria), detritus, and Trichodesmium. In this model, Trichodesmium growth is controlled by the availability of light and its intrinsically slow growth rate. Trichodesmium interacts with no Trichodesmium phytoplankton by competing for light, which depresses its own growth, and by pumping nitrogen into the system, which stimulates no Trichodesmium phytoplankton growth. The ecosystem model has been embedded in a three dimensional circulation model of the Atlantic ocean from 45◦ N to 20◦ S (the Miami Isopycnal Coordinate Ocean or MICOM) (Coles et al., 2004a,b; Hood et al., in press). The coupled physical–biological model reproduces the general characteristics of the observed Trichodesmium biomass distribution. Furthermore, the model generates a clear successional sequence of events including a secondary phytoplankton bloom known as an ‘echo’ bloom. A recent analysis of remote sensing data provides support for the occurrence of such echo blooms fueled by nitrogen fixation (Coles et al., 2004a,b). A recent comparison of the coupled model with and without Trichodesmium suggests that it may be essential to explicitly represent both functional groups, the phytoplankton and the nitrogen fixers, in the ecosystem equations to accurately simulate the spatial and temporal variability of primary productivity and of new nitrogen inputs from nitrogen fixation (Coles et al., 2004a,b). To better understand the biological response of the system to physical forcing, we present in this paper an ecosystem model of reduced dimensionality capturing the essential components of the larger system by Hood et al. (2001). The model consists of three differential equations for the dynamics of phytoplankton, nitrogen fixers, and nutrients. It includes two types of nonlinear interactions between primary producers: the competition of phytoplankton and Trichodesmium for light, and the positive feedback resulting from Trichodesmium recycling.
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
489
Fig. 1. Model structure.
A typical simulation of the model in time, with forcing by a varying mixed-layer depth, reveals a clear successional sequence including a secondary no Trichodesmium phytoplankton bloom. We explain this sequence of events through the stability analysis of the three different steady states of the model. Our analysis shows the existence of a critical biological parameter, the ratio of normalized growth rates, which is non-dimensional and determinines the occurrence of ‘echo’ blooms and the specific sequence of events following a physical perturbation. 2. Model formulation The model consists of three coupled interacting components, describing changes in the concentrations of Nutrient (N) (mmol m−3 ), Phytoplankton (P) (mmol m−3 ) and Trichodesmium (T) (mmol m−3 ). The structure of the model is illustrated in Fig. 1, with arrows indicating flows of matter between the three compartments. The input to the system 0 is D Z (N − N) which represents nutrients entering from below the mixed layer. Arrows not ending in any compartment represent losses from the system. We assume that the mixed layer is thoroughly mixed at all times. Thus, the dynamics of the system can be represented by the three coupled differential equations, dP = U P P − SP P dt dT = G T T − ST T dt D dN = (N 0 − N) − U P P + ηS P P + γ ST T. dt Z We now discuss in detail the different terms and parameters of the model.
(1)
490
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
2.1. Phytoplankton growth The growth rate U P of phytoplankton is a function of light, nutrients and mixed layer depth. If light and nutrient limitation are independent, then we can write U P = I (z, P, T ) × µ P × u(N)
(2)
where I (z, P, T ) is the local light level, µ P is the maximum phytoplankton growth, I is the functional response of growth to light, and u(N) describes nutrient limitation. The local light level varies with z, the depth of the water column and the self-shading by phytoplankton (P), Thrichodesmium (T), and the attenuation by water itself. We then have (Steele, 1962) I = I (z, P, T ) = I0 e−[kx +k P ( P+T )]z
(3)
where I0 is the irradiance intensity at the sea surface, k P is the phytoplankton-specific attenuation coefficient, and k x is the short-wave extinction coefficient. The local light averaged over a mixed layer of depth Z , is given by Z ¯I = 1 I0 e−[kx +k P ( P+T )]z dz (4) Z 0 and I¯ =
I0 1 − e−[kx +k P ( P+T )]Z . Z [k x + k P (P + T )]
(5)
The function u(N) describes the nutrient uptake rate of phytoplankton and satisfies the following general assumptions: it is non-negative and increasing, saturates at high nutrient levels, and vanishes at zero (Hale and Somolinos, 1983). That is, u(N) is a continuous function defined on [0, ∞) satisfying u(0) = 0,
du >0 dN
and
lim u(N) = 1.
N→∞
(6)
In particular, this kind of function includes a Holling’s type II (Holling, 1959) or Michaelis–Menten function (Wroblewski and Richman, 1987) u(N) =
N , N+K
(7)
where K is the half-saturation constant. By neglecting the term e−[kx +k P ( P+T )]Z we write the production of new phytoplankton as N I0 µ P P × . Z [k x + k P (P + T )] N+K
(8)
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
491
2.2. Trichodesmium growth The growth of T is given by I0 µ T T (9) Z [k x + k P (P + T )] which includes the competition with P and other T cells for light through shading. It does not incorporate, however, any nutrient limitation because Trichodesmium fix nitrogen rather than utilize nitrate, and the nitrogen is plentiful in supply. 2.3. Loss and recycling of phytoplankton and Trichodesmium The loss of phytoplankton (Trichodesmium) is modeled by the expression S P P (ST T resp.), where S P and ST are the specific loss rates. These terms represent the mortality of cells due to death and grazing and the loss due to respiration. We assume that only a fraction of the dead phytoplankton (ηS P ) and the dead Trichodesmium (γ ST ) are recycled into dissolved nutrients. The model further assumes that the recycling is instantaneous, thus neglecting the time required to regenerate nutrient from dead biomass by bacterial decomposition. However, a delay in nutrient recycling is present in natural systems and varies inversely with temperature (Whitttaker, 1975). Model (1) can easily be modified to incorporate these delays in the nutrient equation, t D dN = (N 0 − N) − U P P + ηS P f (t − s)P(s)ds dt Z −∞ t g(t − s)T (s)ds. (10) + γ ST −∞
The delay-kernels f (s) and g(s) are non-negative bounded functions defined on [0, ∞) and describe the contribution to the nutrient recycled at time t of phytoplankton and Trichodesmium cells that died in the past. We normalize the kernels such that ∞ ∞ f (s)ds = 1, g(s)ds = 1. (11) 0
0
We note that the presence of the distributed time delays does not affect the equilibrium values. For simplicity, we will assume that the nutrient recycling of the phytoplankton and Trichodesmium is instantaneous, so that f (s) = g(s) = δ(s), the Dirac delta function. 2.4. Mixed-layer depth We assume that the mixed-layer depth Z varies during the year as follows, 0 ≤ t < 60 200 60 ≤ t ≤ 270 Z = 20 2t − 520 270 ≤ t ≤ 360, where t is measured in days and Z , in meters. In the simulations, we consider that the mixed layer-depth varies continuously throughout the year, except for a discontinuous jump of 180 m at t = 60. This piecewise-linear function provides a simple approximation of the three major features of the annual cycle (Evans and Parslow, 1985; Fasham, 1995; Hood
492
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
et al., 2001): a rapid shallowing in the spring, a constant shallow layer throughout the summer, and a deepening during the autumn and winter. 2.5. Mixing The model assumes that the ocean is divided into two completely mixed layers: an upper one, biologically active, which contains nutrients and plankton of constant concentration 0 with depth, and a lower one, inactive, with nutrients but no plankton. The D Z (N − N) term models the exchange of nutrients with the water below the mixed layer, which is assumed to have a constant concentration N 0 . The mixing occurs at a constant rate D (m day−1 ). D Z represents the fraction of the mixed layer that is exchanged with deep water through diffusive processes. 3. Interpretation of the ‘echo effect’ Before proceeding to the mathematical analysis of the model, we illustrate with three representative simulations the biological response of the system to changes of the mixedlayer depth. These simulations correspond to different values of the non-dimensional critical parameter ζ =:
µT ST µP SP
relating the maximum growth rate of T to that of P. The growth rates in ζ are normalized by the respective loss rates of these two groups. We show here that the different responses in the simulations can be explained through the existence of three different steady states whose local stability properties vary with ζ . The existence and stability properties of these steady states are studied in detail in the next section (Section 4). Fig. 2(a) shows a typical simulation of the system for a low growth rate of T relative 0 to that of P (ζ < N 0N+K < 1). For this choice of parameters, we observe the typical successional sequence already described for the higher dimensional model (See Table 1). In this sequence, mixed-layer depth shallowing stimulates a phytoplankton bloom, which is in turn followed by a Trichodesmium bloom after nutrients become depleted. The lack of nutrients and the concomitant decrease of P, releases T from competition allowing it to grow in spite of its maximum low growth rate. The Trichodesmium bloom then stimulates phytoplankton growth and elevates phytoplankton biomass. We refer here to this specific successional pattern as the ‘echo effect’. Fig. 2(b) shows the corresponding phase-space diagram and the three steady states E 0 = (0, 0, N 0 ), E 2 = (P2 , 0, N2 ), and E 3 = (P ∗ , T ∗ , N ∗ ). As the successional pattern proceeds, the trajectory of the system moves from the proximity of E 0 , to that of E 2 , and finally, E 3 . In the winter, the values of P and T remain low and the nutrient concentration approaches that of the sub-mixed-layer water, N 0 . We will show in Section 4 that the steady state E 0 = (0, 0, N 0 ) is globally stable provided the mixed layer is sufficiently deep, as is indeed the case in the winter. In the spring, a bloom of P occurs, as T remains low and N is almost entirely consumed. We will show that, as the mixed layer shallows, the steady state E 2 = (P2 , 0, N2 ) becomes locally stable. For the even shallower summer mixed layer,
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
493
Fig. 2. (a) Temporal dynamics of the model showing an ‘echo’ bloom: phytoplankton (P) (dashed line), Trichodesmium (T) (dashed–dotted line), nutrients (N) (thick line). The mixed-layer depth MLD is in meters (solid line). (b) Phase-space diagram.
E 2 loses stability and E 3 = (P ∗ , T ∗ , N ∗ ) becomes in turn locally stable. The tendency of the system towards E 3 explains the bloom of T and the secondary bloom of P. Finally, by the end of the summer, the deepening of the mixed layer returns the system towards E 0 . Fig. 3 shows the simulated time series for a larger growth rate of Trichodesmium such that N0 < ζ < 1, N0 + K with all other parameters kept at their previous values. In this case, the trajectory in phase space goes from E 0 directly to E 1 = (0, T1 , N1 ), then to E 3 , returning to E 0 . The absence of E 2 which remains unstable at all times, underlies the lack of a spring bloom of P. The system responds to the shallowing of the mixed layer directly with a bloom of T and the echo effect does not occur. There is instead a single P bloom stimulated in part by nitrogen fixation. Fig. 4 shows the simulated times series for an even larger growth rate of Trichodesmium, such that ζ > 1. In this case Trichodesmium blooms in the spring but both phytoplankton blooms are absent. Biologically the phytoplankton spring bloom does not occur because
494
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
Fig. 3. (a) Temporal dynamics with a single phytoplankton bloom following that of the nitrogen fixers. There is no phytoplankton spring bloom. (b) Phase-space diagram (see text).
the Trichodesmium bloom first and then the phytoplankton become light limited, as the result of the negative feedback between these groups. The trajectory in phase space goes from the proximity of E 0 to E 1 = (0, T1 , N1 ), returning to E 0 . The above cases, as well as a more general picture of the system’s dynamics, are summarized in a two-dimensional bifurcation diagram (Fig. 5). This diagram illustrates the stability of the equilibria as a function of two parameters: mixed layer depth, MLD, and the ratio of maximum growth rates ζ . The functions and parameters used to plot the different regions in this graph are given in the next section. To predict the specific response to a shallowing of the mixed layer, a vertical line can be drawn for the relevant value of ζ . The regions crossed by this line tell us whether an echo effect will occur, or alternatively, whether the main bloom of P will be absent. As we will show, however, this graph gives us the global dynamics of E 0 , E 1 , and E 2 but not that of the interior steady state E 3 . To complete the picture, we need to consider the additional parameter µT . Fig. 6 shows that for ζ = 0.1374, the mixed layer fixed at 20 m and all other parameters as in the simulation of Fig. 2 except for an extremely large (and biologically unrealistic) value of µT = 45.9, the trajectory in phase space goes to a periodic orbit shown in Fig. 6. This result, which is relevant to the appearance of multiple blooms, will be proven mathematically at the end of next section.
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
495
Fig. 4. (a) Temporal dynamics with no phytoplankton bloom following the shallowing of the mixed-layer depth. (b) Phase-space diagram.
4. Steady-state analysis System (1) takes the specific form I0 µ P P N dP = − SP P dt Z [k x + k P (P + T )] N + K dT I0 µ T T = − ST T dt Z [k x + k P (P + T )] dN D N I0 µ P P = (N 0 − N) − + ηS P P + γ ST T, dt Z Z [k x + k P (P + T )] N + K
(12)
with initial conditions P(0) = P0 ≥ 0,
T (0) = T0 ≥ 0,
N(0) = N0 ≥ 0.
By adding all the three equations, we can see that d (P + T + N) = −D N − (1 − η)S P P − (1 − γ )ST T dt D I0 µ T T + N0 + Z [k x + k P (P + T )] Z
(13)
496
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
Fig. 5. Bifurcation diagram as a function of mixed-layer depth (MLD) and normalized growth rate ζ (see text).
≤ −D1 (P + T + N) +
I0 µ T D + N0 ZkP Z
where D1 = max D Z , (1 − η)S P , (1 − γ )ST . Thus, the solutions are bounded. From the first equation of system (12), it also follows that P = 0 is an invariant subset, that is, P = 0 if and only if P(t) = 0 for some t. Thus, P(t) > 0. A similar argument holds for T = 0 using the second equation of system (12). We further have N(t) positive. Indeed, if N(t) = 0 for some t, then dN dt > 0 which implies that N remains positive. This proves: Proposition 1. All solutions of system (12) under (13) are bounded and positive. We begin with equilibrium E 0 = (0, 0, N 0 ) which exists for all parameter values. The local stability of E 0 is determined by the eigenvalues of the Jacobian matrix evaluated at this point, I0 µ P N 0 0 0 − SP Z kx N 0 + K I µ 0 T M0 = 0 − ST 0 . Z kx I0 µ P N 0 − γ S −D + ηS P T 0 Z kx N + K The eigenvalues of M0 are I0 µ P N 0 I0 µ T D y2 = − ST , y3 = − . − SP , 0 Z kx N + K Z kx Z Hence, the local stability of E 0 depends on the depth of the mixed layer as described below. This proves: y1 =
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
497
Fig. 6. (a) Temporal dynamics showing periodic solutions. The concentration of T is magnified (×20). (b) Phasespace diagram.
Proposition 2. If I0 µ T Z > max ψ(ζ ), , ST k x
(14)
where ψ(ζ ) =
I0 µ T N0 1 ST k x N 0 + K ζ
then E 0 is asymptotically stable. By changing slightly condition (14) to I0 µ P I0 µ T Z > max , , S P k x ST k x
(15)
we can establish the global stability of E 0 using an a priori estimate of P, T , and N. First we can see that dP I0 µ P P N = − SP P dt Z [k x + k P (P + T )] N + K
498
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
≤
I0 µ P − S P P. Z kx
which using Gronwall’s inequality (Hale, 1969, p. 35), gives I0 µ P P(t) ≤ P0 exp − SP t . Z kx I0 µ P S P kx ,
limt →∞ P(t) = 0. A similar argument follows for T with I0 µ T − ST t . T (t) ≤ T0 exp Z kx
Thus, if Z >
Thus, if Z >
I0 µT S T kx
(16)
, limt →∞ T (t) = 0. On one hand we have
I0 µ P P D N dN = (N 0 − N) − + ηS P P + γ ST T dt Z Z [k x + k P (P + T )] N + K D ≤ (N 0 − N) + ηS P P + γ ST T Z D 0 ≤ (N − N) + (ηS P P0 + γ ST T0 ) exp(−δt) Z
(17)
I0 µ T I0 µ P − SP , − − ST . δ = min − Z kx Z kx
(19)
(18)
where
On the other hand, we also have D N I0 µ P P dN = (N 0 − N) − + ηS P P + γ ST T dt Z Z [k x + k P (P + T )] N + K D I0 µ P ≥ (N 0 − N) − P. Z kx Z
(20)
Let h = N − N 0 . Then, by combining inequalities (18), (20) and using similar arguments to those already used for P, it follows that limt →∞ h(t) = 0. In conclusion, under condition (15), limt →∞ (P(t), T (t), N(t)) = E 0 . This proves: Proposition 3. Under condition (15) E 0 is globally stable. By setting P = 0, ddtT = 0 and for equilibrium E 1 ,
dN dt
I0 µ T = ST Z (k x + k P T ) D 0 (N − N) + γ ST T = 0. Z The solutions are given by T1 =
I0 µ T kx − Z k P ST kP
= 0, we next obtain the following system of equations
(21) (22)
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
N1 = N 0 +
γ ST
D Z
I0 µ T kx − Z k P ST kP
499
.
The steady state N1 is interesting because N1 > N 0 , unlike the case for normal NPZ models (Edwards, 1997; Edwards and Brindley, 1997). Biologically, this difference shows that the presence of Trichodesmium increases the nitrate steady-state value above that for which Trichodesmium is absent, because of nitrogen fixation. The equilibrium E 1 = (0, T1 , N1 ) is positive if and only if Z < kIx0SµTT . Notice that this condition immediately implies the instability of E 0 . Furthermore, as Z approaches kIx0SµTT , E 1 approaches E 0 . The Jacobian matrix at E 1 is µ P ST N1 − SP 0 0 µT N1 + K I0 µ T Z ST k x kPT M1 = − − 1 ST 0 2 Z (k x + k P (P + T )) I0 µ T I0 µ P kx + k P T N D + ηS P − γ ST − Z (k x + k P (P + T ))2 N + K Z with eigenvalues N1 µ P ST y1 = − SP , µT N1 + K
y2 =
Z ST k x − 1 ST , I0 µ T
By writing µT S P , µ P ST
ζ =
it follows that a necessary condition for the stability of E 1 is N0 . N0 + K This proves: ζ >
Proposition 4. If one of the following conditions is satisfied (i)
N0 N 0 +K
< ζ < 1 and ϕ(ζ ) < Z
1
and
Z
N0 , N 0 +K
E 2 is locally asymptotically stable if and only if I0 µ T D ζK − N0 Z< + and ζ > 1. k x ST λk x 1 − ζ
Hereafter we assume that ζ < 1. The interior steady state E 3 = (P ∗ , T ∗ , N ∗ ) is then given by Kζ 1−ζ γ ST I0 µT D 1 P∗ = − k x + (N 0 − N ∗ ) λ + γ ST k P Z ST Z 1 λ I D µ 0 T − k x − (N 0 − N ∗ ) . T∗ = λ + γ ST k P Z ST Z N∗ =
Furthermore the following condition Z < min(ϕ(ζ ), φ(ζ )) ensures the positivity of (P ∗ , T ∗ ), with I0 µ T ζK Dk P − N0 φ(ζ ) = + . ST k x 1−ζ λk x ST The Jacobian matrix at E 3 is given by m 1 m 1 m 13 M = m2 m2 0 m 31 m 32 m 33 where
Z m 1 = −k P S P ST P∗ < 0 I0 µ T Z m 2 = −k P ST2 T∗ < 0 I0 µ T µP K P∗ m 13 = ST >0 µT (K + N ∗ )2
(27) (28) (29)
502
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
m 31 = −(k x + k P T ∗ )ST S P
Z + ηS P I0 µ T
m 32 = −m 1 + γ ST > 0 D m 33 = − − m 13 < 0. Z The corresponding characteristic equation is λ3 + a1 λ2 + a2 λ + a3 = 0
(30)
where a1 = −(m 1 + m 2 + m 33 ) a2 = m 33 (m 1 + m 2 ) − m 31 m 13 a3 = m 2 m 13 (m 31 − m 32 ). By the Routh–Hurwitz criteria, all roots of Eq. (30) have negative real parts if and only if a1 > 0,
a3 > 0
and
a1 a2 − a3 > 0.
We have that a1 > 0 and a3 = m 2 (η − 1) > 0. We also have a1 a2 − a3 = −(m 1 + m 2 + m 33 )(m 33 (m 1 + m 2 ) − m 31 m 13 ) − m 2 (η − 1) = −(m 1 + m 2 + m 33 )m 33 (m 1 + m 2 ) + (m 1 + m 2 + m 33 )m 31 m 13 − m 2 (η − 1) = −(m 1 + m 2 + m 33 )m 33 (m 1 + m 2 ) + (m 1 + m 33 )m 31 m 13 + m 2 (m 31 m 13 − (η − 1)). Thus, if m 31 < 0
and
(m 31 m 13 − (η − 1)) < 0
then a1 a2 − a3 > 0. One can see that if Z > θ (ζ ) then (m 31 m 13 − (η − 1)) < 0, where θ (ζ ) =
I0 µ T k x ST
γ ST (1−η) 2 k P (1 − ζ ) − S P ST (λ + γ ST )K γ ST D 2 0 k P (1 − ζ ) + Z (N + K )
.
The function θ is negative for the model parameters chosen for the main run (Fig. 2). This proves: Proposition 6. The positive steady state E 3 is locally stable if η
I0 µ T < Z < min(ϕ(ζ ), φ(ζ )). k x ST
Now assume that the steady state E 3 is locally stable. We would like to know if E 3 loses its stability when one of the parameters is varied. Note the local stability of E 3 is equivalent to the conditions a1 > 0, a2 > 0 and a1 a2 > a3 . Thus, we would like to know whether
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
503
these conditions will be violated when one of the parameters changes. More specifically, we would like to determine if the real part of the roots of the characteristic equation (30) changes from negative to zero to positive, that is, if a Hopf bifurcation occurs at E 3 . We choose µT , the maximal growth rate of Trichodesmium, as the bifurcation parameter, and address the existence of a critical value µ¯ T such that a1 (µ¯ T )a2 (µ¯ T ) = a3 (µ¯ T ), d [a1 (µT )a2 (µT ) − a3 (µT )] < 0. dµT µT =µ¯ T The characteristic equation (30) becomes (λ + a1 (µ¯ T ))(λ2 + a2 (µ¯ T )) = 0, which has a negative real root λ1 = −a1 (µ¯ T ) and a pair of purely imaginary roots λ2,3 = ±i a2 (µ¯ T ). By continuity, for µT ∈ (µ¯ T − θ, µ¯ T + θ ), with θ > 0 sufficiently small, the eigenvalues are of the form λ1 = −a1 (µT );
λ2,3 = σ (µT ) ± iν(µT ).
To determine if a Hopf bifurcation occurs when µT = µ¯ T , we need to verify the transversality condition. By substituting λ2 into the characteristic equation (30), calculating the derivative with respect to µT , and separating the real and imaginary parts, we obtain dσ (µT ) dν(µT ) − B(µT ) + C(µT ) = 0, dµT dµT dσ (µT ) dν(µT ) B(µT ) + A(µT ) + R(µT ) = 0, dµT dµT A(µT )
where A(µT ) = 3(σ 2 (µT ) − ν 2 (µT )) + 2a1(µT )σ (µT ) + a2 (µT ), B(µT ) = 6σ (µT )ν(µT ) + 2a1 (µT )ν(µT ), C(µT ) = a1 (µT )(σ 2 (µT ) − ν 2 (µT )) + a22 (µT )σ (µT ) + a3 (µT ), R(µT ) = 2a1 (µT )σ (µT )ν(µT ) + a2 (µT )ν(µT ).
Notice that σ (µ¯ T ) = 0, ν(µ¯ T ) = have dRe λ2,3 = σ (µ¯T ) dµ T
µT =µ¯T
= − > 0.
d dµT
√ aµ¯ T . By solving for
dσ dµT
[a1 (µT )a2 (µT ) − a3 (µT )] 4[a12(µT ) + a2 (µT )]
and setting µT = µ¯ T , we
µT =µ¯T
504
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
Table 1 Model parameter symbols, units and values used in the simulations of the echo bloom Description Biological parameters Maximum phytoplankton growth rate Maximum Trichodesmium growth rate Phytoplankton natural mortality rate Trichodesmium natural mortality rate Sat. const. for N uptake by P Phytoplankton-specific attenuation coeff. Recycling efficiency of phytoplankton Recycling efficiency of Trichodesmium Light efficiency at the sea surface Physical parameters Shortwave extinction coefficient Concentration below mixed layer Mixing rate Mixed layer depth
Symbol
Value
Units
µP µT SP ST K kP η γ I0
3.22 0.17 0.08 0.09 0.5 0.0223 0.5 0.8 0.116
day−1 day−1 day−1 day−1 mmol m−3 m2 mmol−1 Dimensionless Dimensionless Dimensionless
kx N0 D Z
0.03 0.6 0.005 20–200
m−1 mmol m−3 day−1 m
Thus, the real part of λ1,2 changes from negative to zero at µT = µ¯ T and becomes positive when µT > µ¯ T . A summary of the above analysis is given in the following proposition on the bifurcation of system (12). Proposition 7. If there exists a positive number µ¯ T such that a1 (µ¯ T )a2 (µ¯ T ) = a3 (µ¯ T )
and
d [a1 (µT )a2 (µT ) − a3 (µT )] µ¯ T , E 3 becomes unstable and a family of periodic solutions bifurcates from E 3 . 5. Discussion We have presented a model with two functional groups, the traditional ‘compartment’ P and an additional one representing nitrogen fixers. The model includes two types of nonlinear interactions between these groups: the competition of phytoplankton and Trichodesmium for light, and the positive feedback resulting from nutrient inputs via Trichodesmium. Although quite simple, this system captures the essential features of the larger ecosystem model introduced by Hood et al. (2001) and used in coupled biological–physical simulations (Hood et al., in press). Through the stability analysis of three different equilibria, we have shown the existence of a critical biological parameter ζ , the ratio of normalized growth rates. A two-dimensional bifurcation diagram based on ζ and mixedlayer depth, was developed to explain the different possible paths of the system in response to the yearly seasonal variation of mixed-layer depth. Mixed-layer depth has long been recognized as a critical physical parameter for the occurrence of phytoplankton blooms. Our
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
505
analysis shows that, when nitrogen fixation is considered, the biological parameter ζ is also important, specifically to the occurrence of a secondary ‘echo’ bloom. As for other models, the first bloom in the spring occurs in response to the physics. The second one is then due to the biological interactions and occurs when ζ is small. It is interesting to note that this second bloom occurs in the autumn. Thus ‘echo’ blooms fueled by nitrogen fixation provide an alternative explanation for phytoplankton increases observed at that time of the year, typically explained by the deepening of the mixed layer while sunlight is still plentiful. A third axis not included in the bifurcation diagram concerns the maximum growth rate of the nitrogen fixers. For fast growth, it is still possible to retain a low ζ provided the mortality rate is also high. Under these conditions, we have shown the tendency of the system to oscillate. Even decaying oscillations would set a scenario in which multiple phytoplankton blooms follow the main spring bloom. Although parameters in this region appear so far unrealistic, the feasibility of this scenario remains open with the recent discovery of small nitrogen fixers in the picoplankton component, whose identity and biology is only beginning to be understood (Zehr et al., 2001). Two potentially important limiting nutrients are not represented in our model: phosphorus and iron. Fe is required for nitrogenase, the enzyme responsible for N2 reduction activity. Because Fe is extremely insoluble in oxygenated seawater, the bioavailability of Fe is controlled by the presence of a variety of organic ligands that enhance Fe solubility (Rue and Bruland, 1995; Karl et al., 2002). Questions therefore arise concerning the proximate control(s) of N2 -fixation rates, including Fe, phosphorus, and other growth factors. Although such controls have been omitted from our equations, our results and general approach may prove useful in interpreting the output of alternative, more complicated models. In particular, it will be interesting to ascertain whether a similar bifurcation diagram underlies the trajectories of alternative models in which a factor other than light controls growth and underlies the competitive interaction of the two functional groups. For environments in which competition for light is not critical (e.g. the Pacific ocean, T. Michaels, personal communication), we propose that the mixed-layer depth axis in the diagram would be replaced by the new limiting factor (such as iron). This possibility remains to be determined. Another extension concerns the treatment of explicit space and the relevant spatial scales at which to address the spatio-temporal dynamics of phytoplankton ecosystems with nitrogen fixation. The large numerical simulations of coupled biological–physical models (Coles et al., 2004a; Hood et al., in press) have not so far resolved space sufficiently to address the effect of eddies, the mesoscale variability or ‘weather’ of the ocean. This question is particularly intriguing given the recent discovery of small nitrogen fixers (Fuhrman and Capone, 2001; Zehr et al., 2001) with potentially faster growth rates, and therefore also faster temporal scales of response to physical forcing. Little is known about their role in the dynamics of plankton ecosystems. Acknowledgements We thank Raleigh Hood for early discussion on the model formulation. This work was partially supported by an NSF Biocomplexity grant.
506
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
References Bates, N.R., 2001. Interannual variability of oceanic CO2 and biogeochemical properties in the Western North Atlantic subtropical gyre. Deep-Sea Res. II 48, 1507–1528. Bates, N.R., Micheals, A.F., Knap, A.H., 1996. Seasonal and interannual variability of oceanic carbon dioxide species at the US JGOFS Bermuda Atlantic Time-series Study (BATS) site. Deep-Sea Res. 43, 347–383. Bates, N.R., Takahashi, T., Chipman, D.W., Knap, A.H., 1998. Variability of pCO2 on diel to seasonal timescales in the Sargasso Sea near Bermuda. J. Geophys. Res. 96, 15567–15585. Capone, D.G., 2000. The marine nitrogen cycle. In: Kirchman, D. (Ed.), Marine Microbial Ecology. John Wiley, NY, pp. 455–493. Capone, D.G., Carpenter, 1982. Nitrogen fixation in the marine environment. Science 217, 1140–1142. Capone, D.G., Subramaniam, A., Montoya, J.P., Voss, M., Humborg, C., Johansen, A.M., Siefert, R.L., Carpenter, E.J., 1998. An extensive bloom of the N2-fixing cyanobacterium, Trichodesmium erythraem in the central Arabian Sea. Mar. Ecol. Prog. Ser. 172, 281–292. Capone, D.G., Zehr, J.P., Pearl, H.W., Bergman, B., Carpenter, E.J., 1997. Trichodesmium, a globally significant marine cyanobacterium. Science 276, 1221–1229. Coles, V.J., Hood, R.R., Pascual, M., Capone, D.G., 2004a. Modeling the impact of Trichodesmium and nitrogen fixation in the Atlantic ocean. Ocean. J. Geophys. Res. 109, C06007. Coles, V.J., Wilson, C., Hood, R.R., 2004b. Remote sensing of new production fuelled by nitrogen fixation. Geophys. Res. Lett. 31, L06301. Edwards, A.M., 1997. A rational dynamical-systems approach to plankton population modeling. Ph.D. Thesis, University of Leeds. Edwards, A.M., Brindley, J., 1997. Oscillatory behavior in a three-component plankton population model. Dyn. Stability System 11, 347–370. Evans, G., Parslow, J., 1985. A model of annual plankton cycles. Biol. Oceanogr. 3, 327–347. Fasham, M.J.R., 1995. Variations in the seasonal cycle of biological production in subarctic oceans: a model sensitivity analysis. Deep-Sea Res. 42, 1111–1149. Fuhrman, J.A., Capone, D.G., 2001. Nifty nanoplankton. Nature 412, 593–594. Gruber, N., Sarmiento, J.L., 1997. Global pattern in marine nitrogen fixation and denitrification. Global Biogeochem. Cycles 11, 235–266. Hale, J.K., 1969. Ordinary Differential Equations. Wiley-interscience, New York. Hale, J.K., Somolinos, A.S., 1983. Competition for fluctuating nutrients. J. Math. Biol. 18, 255–280. Holling, C.S., 1959. Some characteristics of simple types of predation and parasitism. Canad. Entomol. 91, 385–398. Hood, R.R., Capone, D.G., Olson, D.B., 2001. Modeling seasonal and interanual biogeochemical and N2 cycles at BATS. Deep-Sea Res. II 48, 1609–1648. Hood, R.R., Coles, V.J., Capone, D.G., 2004. Modeling the Distribution of Trichodesmium and Nitrogen Fixation in the Atlantic, Ocean. J. Geophys. Res. (in press). Karl, D. et al., 2002. Dinitrogen fixation in the world’s oceans. Biogeochemistry 57, 517–519. Karl, D., Letelier, R., Hebel, D.V., Bird, D.F., Winn, C.D., 1992. Trichodesmium blooms and new nitrogen in the North Pacific gyre. In: Capone, D.G., Reuter, J.G. (Eds.), Marine Pelagic Cyanobacteria: Trichodesmium and other Diazotrophs. Kluwer Academic Publishers, Boston, pp. 219–237. Michaels, A.F., Olson, D., Sarmiento, J., Ammerman, J., Fanning, K., Jahnke, R., Knap, A.H., Lipschultz, F., Prospero, J., 1996. Inputs, losses and transformations of nitrogen and phosphorus in the Pelagic North Atlantic Ocean. Biogeochemistry 35, 181–226. Rue, E., Bruland, K., 1995. Complexation of iron (III) by natural organic ligands in the Central North Pacific as determineted by new competitive ligand equilibration/adsorptive cathodic voltametric method. Mar. Chem. 117–138. Sambrotto, R.N., Savidge, G., Robinson, C., Boyd, P., Takahshi, T., Karl, D.M., Langdon, C., Chipman, D., Marra, J., Codispoti, L., 1993. Elevated consumption of carbon relative to nitrogen in the surface ocean. Nature 363, 248–250. Steele, J.H., 1962. Environmental control of photosynthesis in the sea. Limnol. Oceanogr. 7, 137–150. Whitttaker, R.H., 1975. Communities and Ecosystems. Macmillan, New York.
K. Boushaba, M. Pascual / Bulletin of Mathematical Biology 67 (2005) 487–507
507
Wroblewski, J.S., Richman, J.G., 1987. Nonlinear response of plankton to wind mixing events—implications for survival of larval northern anchovy. J. Plank. Res. 9, 103–123. Zehr, J.P., Waterbury, J.B., Turner, P.J., Montoya, J.P., Omoregie, E., Steward, G.F., Hansen, A., Karl, D.M., 2001. Unicellular cyanobacteria fix N2 in the subtropical North Pacific Ocean. Nature 412, 635–638.
