- Email: [email protected]
A coral bleaching model
Nonlinear Analysis: Real World Applications 16 (2014) 65–73
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
A coral bleaching model Peter L. Antonelli a,∗,1 , Solange F. Rutz b , Paul W. Sammarco c , Kevin B. Strychar d a
Department of Mathematical Sciences, University of Alberta, Edmonton, AB, Canada
b
Federal University of Pernambuco - UFPE, Mathematics Department, Recife, PE, Brazil
c
Louisiana Universities Marine Consortium - LUMCON, Chauvin, LA, USA d Annis Water Resources Institute, Grand Valley State University, Muskegon, MI, USA
article
info
Article history: Received 25 May 2013 Accepted 14 September 2013
abstract Volterra–Hamilton systems theory is used to model coral bleaching. The algal and coral partners in the obligate symbiont organism exchange compounds, each one producing compounds benefiting the other. It is found that the production equations, based on a homogeneous cost functional, have unstable solutions which become stable in a constant temperature environment. Temperature spikes and global warming occur in the noisy background. The variety of clade distributions for the algal symbiont is shown to be a consequence of Nelson’s stochastic mechanics. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Coral bleaching is a result of the expulsion of symbiont algae, known as zooxanthellae, by their host organisms commonly called stony corals, or hermatypic scleractinians, the main reef builders [1]. Bleaching is a major problem globally because it threatens the existence of tropical reefs [2–4]. This process is called bleaching because of the loss of pigmentation contributed to the holobiont by the zooxanthellae, resulting in the white calcium carbonate exoskeleton of the coral host showing through the still intact but transparent coral tissue. Zooxanthellae exist in a large number of clades or variations of this alga and are found in the interstitial spaces of the coral tissue in all hermatypic or reef-building corals [5,6]. They are also found in soft corals, or octocorallians [7]. Zooxanthellar clades vary between coral species. There is always one clade of zooxanthellae in a given coral species in a given location and time and sometimes more than one clade. Corals can only survive and grow within a relatively narrow range of seawater temperatures and temperature increases predicted to occur in the world’s oceans could conceivably drive today’s reefs towards extinction [2]. This is the ecological problem addressed in this paper. Recent research has shown that the scleractinian host corals have adapted, over evolutionary time, to moderate fluctuations in seawater temperature, exhibiting a relatively broad range of temperature tolerances, while their symbiotic zooxanthellar partners have not [8]. Indeed, the zooxanthellae seem to be catching up, adapting more rapidly to temperature-based environmental fluctuations, which are increasing with time because of global warming [9–11]. The generation time (i.e., length of time it takes to reach age of first reproduction) for corals is on the order of years, while generation time for the zooxanthellae is on the order of days to weeks [11]. This enables the zooxanthellae to reproduce at a considerably greater rate than their coral hosts by possessing a much higher frequency of reproduction. A shorter generation
∗
Corresponding author. Tel.: +1 780 438 7886. E-mail addresses: [email protected] (P.L. Antonelli), [email protected] (S.F. Rutz), [email protected] (P.W. Sammarco), [email protected] (K.B. Strychar). 1 Visiting Professor at UFPE, Recife, PE, Brazil. 1468-1218/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.nonrwa.2013.09.006
66
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
time implies a better ability of these types of populations to adapt, merely on the basis of frequency of reproduction, all other things being equal. But, when generation time is compared to mutation rate in these organisms, it was determined that, on a per generation basis, zooxanthellae actually have a mutation rate lower than expected for an organism of their size [11]. This is counterintuitive, since mutations are one of the factors that drive adaptation and evolution. If this is the case, the question is raised, which of the two organisms comprising the holobiont, is actually responsible for most of the adapting? Could it be that the zooxanthellar mutation rate is sufficient to produce, on an evolutionary time scale, the many types of zooxanthellar clades now known to exist? We modelled these traits in order to understand and possibly predict survivorship and adaptation in corals under environmental conditions of rising seawater temperatures, known to be occurring at present. Our approach is both ecological and physiological in perspective. 2. Basics of our approach There are tens of algal clades which make up zooxanthellae, but for a coral species at a given reef location and time, only a few algal types at most are housed, and these can differ genetically and often in pigmentation. It has been suggested that the relatively short generation time of the zooxanthellae, together with natural selection, would provide a great number of variable zooxanthellar types, on an evolutionary time scale. Under normal environmental conditions, there would be a well-defined average value for their spatial density in the coral tissue, with a relatively small variance over time. We denote this density as, N 1 (t ), with t in a fixed interval [t0 , t1 ], in a single coral colony at a fixed location. Although there may be more than one type of zooxanthella living in symbiosis within a single coral, for simplicity in model building, we will assume here that there is just one. There is considerable variety and turnover of zooxanthellae in a coral—based upon repeated colonization. We presume this colonization process is controlled by interactions between the coral cells and the symbionts via an immune response, which is in turn genetically determined but may be influenced by the local, proximal, ecological and environmental conditions. This interaction can be described as in classical mathematical ecology using approximately constant coefficients in first order nonlinear differential equations [12–15]. There are two populations, one of coral cells, the other of algal cells, and these are taken to be weakly coupled which is expressed by coefficients being constant, or nearly so. This is in contrast to the strongly coupled case in which the coefficients actually depend explicitly on the x1 , x2 (amounts produced by, and exchanged with, each partner in a holobiont organism) and even on the ratio of the N [16,17]. Strong couplings do occur in nature and involve robustly stable, genetically programmed, chemical exchanges [18–20]. The term ‘‘robustly stable’’, as used here, means that the production stability measure, K , is positive and stays positive both before and after small perturbations of the x-values are performed in the expression for K . However, the bleaching situation addressed here seems not to be robustly stable, so there must be some small perturbations of x-values which render K non-positive. A specific zooxanthellar clade in residence in a specific coral tissue is not considered a permanent resident, despite the obligate relationship between the host and the symbiont. It might be referred to as a ‘‘fair weather symbiont’’ and a clade within a coral may change over time. This is far different from the obligate symbionts in, say, root fungi in vascular plants, or organelles, as mitochondria and chloroplasts in eukaryotes [18,19]. The coral acquire zooxanthellae as symbionts and these microalgae enable the corals to produce aragonitic calcium carbonate, the basic substrate of a coral reef. This ability is an evolutionary innovation which came about during the late Triassic period [21,22]. Before this, corals had no skeletons and did not form reefs (e.g. rugose corals). Nowadays, of course, each population in the holobiont gains from the association and some in ways naturally described as allometric. For example, the coral host #1, obtains photosynthates (i.e. sugars and oxygen) from the algae in an amount denoted x1 = b1 ln(m1 ) + c1 , while the algae, #2, receive nitrogenous wastes (e.g., CO2 and water) from the coral in the amount, x2 = b2 ln(m2 )+ c2 . Here, the b and the c are allometric constants and we say that the x are surrogates of biomass. Our model applies to all surrogates of biomass, not just this particular interpretation involving sugars and nitrogen. Further, we consider this obligate relationship as relatively weak and not linked or locked. Zooxanthellar turnover does certainly occur in their natural habitat. Indeed, corals expel zoox, mostly already dead—those which were not able to survive the rise in sea water temperature. One outstanding question is, does the holobiont organism replace expelled algae with a new viable zooxanthellar clade in time to permit continued reef calcium carbonate production? We believe that it does because it follows from the mathematical model below. We will assume holobiont production makes efficient use of available energy (ATP) for the system as a whole; so we will adjoin Volterra production equations, dxi /dt = k(i) N i to the ecological equation system of first order non-linear differential equations, which are to be chosen below [15,17]. This is the basic idea of analytical trophodynamics theory, which focuses on stability properties of the production processes [23,17]. Again, population #1 is the zooxanthellae and #2 is the coral species polyp population. (The parentheses on the lower index indicate that the equation is not summed from 1 to 2. However, the Einstein summation convention, summing on repeated upper and lower indices, is used throughout.) The parameters ki are positive proportionality constants and can be set equal to unity for convenience in this particular model. Since the combined system of equations is of second order, we can use the calculus of variations and implement simple regularity conditions which ensure physiological production processes
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
67
minimize production cost (ATP). We choose a cost function F which is positive and allow for explicit dependence on the produced substances, x1 , x2 , their production rates dx1 /dt, dx2 /dt, and on clock time, t. Moreover, the optimal growth of the two partner populations would, under ideal conditions, be in synchrony – strongly correlated – in order to ensure the spatial distribution of zooxanthellae within the coral cell population remains approximately constant. In reality, this may not be the case because of external factors which can alter the zooxanthellar populations (e.g., disease, ocean acidification, etc.). Here, we will assume that both N 1 and N 2 , are increasing initially at a fixed mean rate, λ. This rate is not to be confused with generation time. Moreover, the measure of joint production will be the cost itself. Thus, the infinitesimal relation ds = F (x1 , x2 , dx1 , dx2 , t , dt ) implies that small positive increments dx1 and dx2 , occurring in time interval, [t , t + dt ], will result in small positive increments ds, in joint production. We will further suppose that twice (or any positive multiple, C of) an increment, would multiply by a factor of two (or by a factor, C ), the joint production in that interval. This is known as the homogeneity assumption, and is biologically sound and has been used previously [24]. It should be noted that the parameter s is akin to size. It is independent of units used to describe it. Indeed, s is to be thought of as a geometric quantity. Thus, a smooth curve γ in the production space has a finite length s given by the integral of ds along γ from one point to another. The positive quantity s is the measure of total joint production along γ , for the entire holobiont population under study. From calculus of variations, the path integral of ds will depend on the smooth path chosen. Generally, we will consider those with dx1 /dt , dx2 /dt both positive and differentiate the integral of F with respect to the path and set it equal to zero. This preliminary will ensure efficiency of the physiologically based production process, because the usual technical assumption of regularity results in equations for the path of minimum cost, between (xio , to ) and (xi1 , t1 ), relative to the given cost functional F . These are the so-called Euler–Lagrange equations associated with F . However, as in the above discussion, these equations must have constant coefficients or nearly so, and this requirement strongly constrains the possible choices of F . There are just three possibilities according to the classification theorem of Antonelli and Matsumoto (see appendix, [25]). However, two of these cost functions are not appropriate, because, in one case, it leads to asymmetrical and biased interactions, like parasitic or commensal ones, and in the other, there are no interactions at all, only two independent logistic equations, one for each population. Only one of the three cost functions leads to an a priori unbiased treatment description of the two cell populations. Namely, ds = [(dx1 )2 + (dx2 )2 ]1/2 exp{(L2 + 1)αi xi + λt /2 + L arctan(dx1 /dx2 ) + ψ(x)}, with L as defined below and ψ(x) = (1/2)[ν1 .(x1 )2 + ν2 (x2 )2 ] [25]. Here, all coefficients denoted by the α , the ν and λ are positive constants. The corresponding Euler–Lagrange equations split naturally into Volterra’s production equations above, with ki = 1, ψ(x) = 0 and ecological equations dN 1 /dt = λN 1 − 2(α2 )N 1 N 2 − (α1 )[(N 1 )2 − (N 2 )2 ], dN 2 /dt = λN 2 − 2(α1 )N 1 N 2 − (α2 )[(N 2 )2 − (N 1 )2 ]. Of course, substituting the Volterra production equations back into these, results in the geodesic equations of the Finsler geometry automatically determined by the cost F . However, the path parameter here, is clock time t, not parameter s (see [15], p. 47). The α , along with λ, determine the conceptually deeper (geometric) structures of the cost F . Note that if N 2 is set to zero in equation #1, then a logistic differential equation is obtained for N 1 , so the algal population will decrease to a (relatively) small size (we assume small carrying capacity K1 ) and will die out from noisy environmental fluctuations alone. Likewise, a logistic is obtained from equation #2, if N 1 = 0 and N 2 will decrease to (relatively) small size (we assume a small carrying capacity K2 ) and will be wiped out by noisy fluctuations. Note that we take, α1 = λ/K1 and α2 = λ/K2 , where K1 and K2 are the (relatively) small carrying capacities for independently living populations of algal and coral cells, respectively. This smallness expresses their relative inability to survive on their own. But, they are partners in a symbiosis that is obligate and this ensures their survival. Remark. Do not confuse ki with Ki or K . They are very different! There is a unique steady-state in the positive (N 1 , N 2 ) quadrant when ψ(x) = 0. This is easy to compute and we leave it for the reader. The constant L = (α2 − α1 )/(α1 + α2 ) may be zero, negative or positive, but the production process is always unstable, since ψ(x) = 0 implies, K = 0, from the K -formula given below. Therefore, it follows there are paths in production space, close to start with, but which diverge eventually. Positive curvature, on the other hand, ensures that all initially close paths remain close. Negative or zero curvature forces instability so that at least some holobiont populations have divergent production curves. If K is strictly less than zero throughout production space, all paths must eventually diverge. Divergent paths, however, are unreliable from an evolutionary viewpoint since reproduction would be unstable, it being just a special case of production in our model, since the x-variables are surrogates of biomass with an assumed allometric relationship between reproductive and total biomass. Remark. It is of interest to write the parameters L and L2 + 1 in terms of the carrying capacities K1 and K2 . They are, L = (K1 − K2 )/(K1 + K2 ) and L2 + 1 = 2(K12 + K22 )/(K1 + K2 )2 . It is obvious that if the zooxanthellae carrying capacity K1 is less than K2 of the coral, then L is negative in the expression for K , while the reverse is true when K2 < K1 . Some standard references for Finsler geometry are [26,17,15,27,25].
68
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
3. Motivation for the model For the unstable case, we have ψ(x) = 0. Yet, certain negative constant terms representing temperature, added to the right-hand sides of the ecological equations, convert the production process to one that is stable. We prove this below using KCC-theory [20]. But, before we do this, we clarify several points: Firstly, if we had positive curvature K (and therefore production stability), then the holobiont would be highly integrated, physiologically. The reason is that positivity of curvature K can only be achieved by including ψ(x) = (1/2).[ν1 (x1 )2 + ν2 (x2 )2 ] in the exponential function and, as a consequence, the coefficients in the ecological equations must explicitly depend on xi , and the interactions, consequently, would exhibit strong coupling. Although there is solid evidence that a chemical exchange between the two components of this symbiotic system is part of an evolutionary innovation for the holobiont, the chemical exchanges occurring are weak compared to the highly integrated robustly stable exchanges mentioned above. Nevertheless, they are obligate, because the hermatypic corals cannot survive without their zooxanthellae, and conversely. Secondly, because there are several clades of zooxanthellae available for the replacement of one currently held by a given species of coral, when and if a seawater temperature spike occurs, and since the (joint) production in the holobiont is considered weak, any stability must be due of the influence of the ambient environment. This is demonstrated below. Third and finally, the model constructed here predicts that (joint) production stability in the holobiont can be achieved, but only over evolutionary time-scales and resulting from physiological changes and integration of the resulting chemical processes into the holobiont genome. Furthermore, the mathematical model predicts that the environmentally forced turnover of zooxanthellar types within a single intact coral colony, does take place, and is caused by rising seawater temperatures and temperature spikes. This is a conclusion from the mathematical model resulting from environmental noise having been added to the model according to the stochastic theory for conservative potential driven systems [28,29]. 4. Constant temperature case The stability measure K is actually the Gauss–Berwald curvature K of the ‘‘production space’’ for the holobiont. The stability scalar K has been computed, via S.F. Rutz’s software package FINSLER [30] based on MAPLE [31]. She expresses it as K = −2(ν1 + ν2 )(L2 + 1) exp{[2φ(x) + 2L arctan(N 1 /N 2 )]},
φ(x) = ψ(x) + 2(L2 + 1)[αi xi ], L = (α2 − α1 )/(α1 + α2 ), L2 + 1 = [(α1 )2 + (α2 )2 ]/(α1 + α2 )2 . It is significant that L and K together are sufficient to uniquely specify the Finsler geometry of the production surface [26,17]. Let us discuss the N-ratio effect on the curvature K . If ν1 and ν2 were positive and, if in addition, N 1 /N 2 was small, it follows that arctan(N 1 /N 2 ) is positive and close to zero which then implies that K , being negative, would be closer to zero. That is, the process of production, although unstable, progresses towards stability. Thus, the reapportionment of the total energy available to the holobiont to include relative population size, indicated by N 1 /N 2 , would help to stabilize production. The importance of this for holobiont evolution is clear since stable production is a necessary condition for stable reproduction. Let us explain further. It is standard procedure in Volterra–Hamilton systems theory (i.e., analytical trophodynamics) to assume that x1 , x2 are surrogates of total accumulated biomass [23,17]. This allows allometric quantification of the products, for example, seed production in plants, fertilized egg production in animals, production of chemicals during growth or even morphological characters. Thus, reproductive materials would be in allometric proportion to the total biomass produced. This is familiar to biologists in the guise of linear regression in log–log plots of phenotypic characters, going back to J. Needham and J. Huxley in the 1930s [32,33]. We argued above that incorporation of negative quadratic terms in φ(x), which necessarily resulted in stability of production, is not well founded. On similar grounds, inclusion of positive versions of φ(x) is not well founded, either. Therefore, we must take both ν1 and ν2 to be negligibly small or zero. This brings us right back to the system with the α restricted to be more or less constant, and for which small perturbations in x-values may alter the stability. 5. KCC-theory proof of environmentally induced stability Let (x1 , . . . , xn ) = (x), (dx1 /dt , . . . , dxn /dt ) = (dx/dt ) = (˙x), and t be 2n + 1 coordinates in an open connected subset Ω of the Euclidean (2n + 1)-dimensional space Rn × Rn × R1 . Suppose that we have d2 xi dt 2
+ g i (x, x˙ , t ) = 0,
i = 1, . . . , n
(1)
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
69
for which each g i is C oo in a neighbourhood of initial conditions ((x)0 , (˙x)0 , t0 ) ∈ Ω . The intrinsic geometric properties of (1) under non-singular transformations of the type
x¯ i = f i (x1 , . . . , xn ), t¯ = t ,
i = 1 , . . . , n,
(2)
are given by the five KCC-differential invariants, named after D. Kosambi [34], E. Cartan [35], and S.S. Chern [36], given below. Let us first define the KCC-covariant differential of a contravariant vector field ξ i (x) on Ω by
Dξ i
=
dt
dξ i dt
1
+ g;ir ξ r ,
(3)
2
where the semi-colon indicates partial differentiation with respect to x˙ r , and we made use of the Einstein summation convention on repeated indices. Using (3), Eq. (1) becomes
D x˙ i dt
= ϵi =
1 2
g;ir x˙ r − g i ,
(4)
defining the first KCC-invariant of (1), the contravariant vector field on Ω , ϵ i , which represents an ‘external force’. Varying trajectories xi (t ) of (1) into nearby ones according to x¯ i (t ) = xi (t ) + ξ i (t )η,
(5)
where η denotes a parameter, with |η| small and ξ (t ) the components of some contravariant vector defined along x = x (t ), we get, substituting (5) into (1) and taking the limit as η → 0, i
d2 ξ i dt 2
+ g;ir
dξ r dt
i
+ g,ir ξ r = 0,
i
(6)
where the comma indicates partial differentiation with respect to xr . Using the KCC-covariant differentiation (3) we can re-express this as
D 2ξ i dt 2
= Pri ξ r ,
(7)
where
Pji = −g,ij −
1 2
g r g;ir ;j +
1 2
x˙ r g,ir ;j +
1 4
g;ir g;rj +
1 ∂ 2 ∂t
g;ij .
(8)
The tensor Pji is the second KCC-invariant of (1). The third, fourth and fifth invariants are:
1 i i i Rjk = (Pj;k − Pk;j ) , 3 i Bjkl = Rjki ;l , i Djkl = g;ij;k;l .
(9)
The main result of KCC-theory is the following: Two systems of the form (1) on Ω are equivalent relative to (2) if and only if the five KCC-invariants are equivalent. In particular, there exists coordinates (¯x) for which the g i (¯x, x˙¯ , t ) all vanish if and only if all KCC-invariants are zero. The tensor D vanishes if and only if g i is quadratic in x˙ , in the case when the first KCC-invariant vanishes. Let us now introduce the notion of a n-dimensional Finsler space as a manifold where, given a coordinate system (x) and a curve xi = xi (t ), the norm of a tangent vector x˙ i to the curve at each point P on xi (t ) is given by the positive metric function F , |˙xi | = F (x, x˙ ), where F is positively homogeneous of degree 1 in x˙ i . From F , a metric tensor is defined as gij (x, x˙ ) = (∂ 2 F 2 /∂ x˙ i ∂ x˙ j )/2, which must be regular in an open region of the tangent bundle, the collection of all tangent vectors to the manifold, and which excludes the origins. The use of the calculus of variations for F leads to (1) with g i (x, x˙ , t ) = γjki (x, x˙ ) x˙ j x˙ k , where the γjki are the Levi-Cività symbols for the Finsler metric tensor gij (x, x˙ ). Berwald’s Gaussian curvature K for two-dimensional Finsler spaces is defined from his famous formula [17] i Rjk = F K mi (lj mk − lk mj ),
(10)
i where Rjk is given by the first equation in (9), li = x˙ i /F is the unit vector in the x˙ i direction, and mi the unique (up
to orientation) unit vector perpendicular to li . Lowering the index on mi via the metric tensor gives mi , which satisfies F (x, m) = gij (x, x˙ ) mi mj ≡ mi mi = 1. If our curvature K is bigger than zero everywhere, then trajectories oscillate back and forth, crossing the reference trajectory. In this case, we say (1) is Jacobi stable. If K ≤ 0 everywhere, trajectories diverge and system (1) is Jacobi unstable [26,15]. This notion of stability is a Lyapunov notion, but it is a whole trajectory concept.
70
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
We consider the RHS of the production equations above with L not necessarily zero. The RHS will be of the form,
−∇ i (V (x)), where V depends only on x-values and the gradient is of the form: g ij (x, dx/dt )∂j with gij = 1/2(∂˙i )(∂˙j )F 2 . The resulting equations are Euler–Lagrange equations for the Lagrangian
L = (1/2)gij (dxi /dt )(dxj /dt ) − V (x). The so-called metric tensor gij depends on the ratios N 1 /N 2 = dx1 /dx2 and is globally defined on the symmetric tensor product T (M ) ⊗ T (M ) with T (M ) denoting the tangent bundle of positive conical regions of the tangent planes, all with origin removed. This is why the non-zero rate condition must be implemented (see [26] vol. 1, p. 570). The deviation curvature tensor for the conservative dynamics of L is found to be
−Pji = Rjki .yk + ∇ i (Vj (x)), V (x) = a/2 exp[2αi xi ],
Vj (x) = ∂j V (x),
g ki = exp[−2αm xm ]h(y)ki
∇ k (Vj (x)) = a.h(y)ki [2(αi )(αj ) − (αr )Grij ] := a.hki Qij . Here, hki depends only on yj = dxj /ds and is positive definite on the tangent cones, and furthermore, Giii = αi ; Giij = αj , if i ̸= j; Gijj = −αi , if i ̸= j. These three-index quantities are symmetric in the two lower indices. It is seen that for all four choices of i, j Qij = [(α1 )2 + (α2 )2 ], so that the gradient of V (x) is positive definite. Of course, Rijk is vanishing as follows from the formula relating it to K . We can now see that the deviation vector ξ i , which is orthogonal to dxi /ds of γ , is oscillating back and forth across the reference trajectory γ in production space and therefore the production process in terms of the allometric x-variables is stable. 6. Nelson’s stochastic mechanics and zooxanthellar clades Define on the Euclidean 2-plane, R2 , the potential V = (a/2) exp[αi xi ]. In fact,
−a(α i ) = −g ij (x)(∂j )(V (x)) which are the right-hand sides of the Euler–Lagrange equations for Lagrangian L = (1/2)gij yi yj − V (x), where yi = dxi /ds and g ij = exp{−2αk xk }. δ ij is the inverse of gij , a Riemannian metric tensor. The Euler–Lagrange equations are dx1 /ds = y1 ,
dy1 /ds = −2(α2 )y1 y2 + (α1 )[(y2 )2 − (y1 )2 ] − a(α1 ),
dx2 /ds = y2 ,
dy2 /ds = −2(α1 )y1 y2 + (α2 )[(y1 )2 − (y2 )2 ] − a(α2 ).
The stochastic theory used here was developed by Edward Nelson of Princeton University and his students for a stochastic approach to derivation of Schrödinger’s equation in classical quantum mechanics. The paths are nowhere differentiable but are continuous and Markov. This theory is mathematically difficult, but has been much studied by mathematical physicists [29]. However, a surprising biological application to social behaviour was presented in [28]. Here, we show how our bleaching model can be transformed into the double harmonic oscillator considered in [28]. A discrete spectrum eigenvalue problem, with known solutions, is obtained. Transforming coordinates, spatial distributions for successively higher total energy are obtained with each a sort of ‘‘mountain range’’ with the lowest energy level yielding a single Gaussian distribution in the phenotypic character space spanned by our surrogates of biomass variables (see appendix, [15]. There are nodal lines defined by setting coordinates equal to zero, and the enclosed regions are isolated from each other in the sense that crossing a nodal boundary, by any sample path starting inside, has probability zero. However, the processes are ergodic within these compartments [15,29,28]. Theorem. There is a non-singular coordinate transformation from x1 , x2 to polar coordinates r , θ given by r (x1 , x2 ) = exp(α1 )x1 + (α2 )x2 ,
θ (x1 , x2 ) = (α2 )x1 − (α1 )x2 , so that our environmentally perturbed system becomes d2 r (s)/ds2 − r [dθ (s)/ds]2 = −kr (s), d2 θ (s)/ds2 + (2/r (s))[dr (s)/ds][dθ (s)/ds] = 0, where k = a[(a1 )2 + (a2 )2 ] and s(t ) = A − B exp{−(λ)t } are convenient definitions of constants and path parameter s in terms of clock time t. This is the polar form of the double harmonic oscillator.
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
71
Fig. 1. Density ρ0 (x, y).
The Newton–Nelson formula [26, page 269] is equivalent to Schrödinger’s equations for the two dimensional harmonic oscillator, but Planck’s constant is replaced by ν 2 , and the stationary solution ψ satisfies
−
ν2 2
∆+V
ψ = λψ,
√ ρ . According to well-known sources, the discrete eigenvalues satisfy √ λn = ν 2 a (n + 1), n = 1, 2, . . . .
where ψ =
The stationary solution for the smallest eigenvalue λ0 is
λ0 = C0 exp −
1 2
a
ν2
r
2
.
The corresponding density function is a single ‘‘bump’’. For the second eigenvalue λ1 there are two independent stationary solutions (see Figs. 1 and 2 and [37] for more graphs), and for the third, λ2 , there are three, and so forth. Finally, the densities, as pictured, are given by squaring the stationary solutions. These densities describe chemical morphology of the holobiont. The upshot of this method of noise modelling is that the total energy (ATP plus increments produced by the gradient of an environmental potential V ) yields distributions (the diagrams above) for the holobiont population, depending on the total energy level, with the larger energies yielding more peaks (more ‘‘mountains in the mountain range’’). Most of the holobiont variety, as discussed above, is understood to be the various clades of the symbiont algae. Furthermore, when a random (seawater temperature) spike event occurs in the environment of a holobiont state situated near a compartmental boundary, the occupied state can continuously move towards the nearby boundary and pass through to a new distribution. Yet, all compartments at a given total energy level are isolated from each other in Nelson stochastic theory. By allowing higher and higher total energy levels, all holobiont types occurring in nature can, in principle, be represented. A random seawater temperature spike can force communication between the isolated compartments. As time goes on, more and more ‘‘mountains’’ are occupied providing the necessary rich assortment of types. To be precise, the parameter a is assumed to increase or decrease by jumps of positive or negative magnitude δ , at random times, the number J+ (S ) of jumps by +δ and the number J− (S ) of jumps by −δ up to time s, being Poisson processes independent of each other and of the background noise in the environment. The requisite computer simulations were published previously [37]. The standard approach to adding noise to Finslerian production equations has been developed for several distinct two-dimensional geodesic systems [25]. However, the Nelson theory for conservative systems has not been developed and not for lack of trying! Fortunately, we are able in this present model to use the Riemannian version successfully. Detailed proofs have appeared in several places, in particular, [15, appendix], [25,26], or [37].
72
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
Fig. 2. Density ρ01 (x, y).
Acknowledgement FACEPE grant BFP-0043-1.01/11. References [1] C.E. Birkeland, The Life and Death of Coral Reefs, Springer, New York, 1997, p. 536. [2] O. Hoegh-Guldberg, Climate change, coral bleaching, and the future of the world’s coral reefs, Mar. Freshw. Res. 50 (1999) 839–866. [3] M.J.H. Van Oppen, J.C. Mieog, C.A. Sanchez, K.E. Fabricius, Diversity of algal endosymbionts (zooxanthellae) in octocorals: the roles of geography and host relationships, Molec. Ecol. 14 (2005) 2403–2417. http://dx.doi.org/10.1111/j.1365-294X.2005.02545.x. [4] M.J.H. Van Oppen, J.M. Lough (Eds.), Coral Bleaching: Patterns, Processes, Causes, and Consequences, in: Ecol. Studies, vol. 205, Springer, New York, 2009, p. 178. [5] R. Rowan, Review—diversity and ecology of zooxanthellae on coral reefs, J. Phycol. 34 (1998) 407–417. http://dx.doi.org/10.1046/j.15298817.1998.340407.x. [6] A.C. Baker, Flexibility and specificity in coral-algal symbiosis: diversity, ecology, and biogeography of symbiodinium, Ann. Rev. Ecol. Evol. Systematics 34 (2003) 661–689. http://www.jstor.org/stable/30033790. [7] O. Hoegh-Guldberg, P.J. Mumby, A.J. Hooten, R.S. Steneck, P. Greenfield, E. Gomez, C.D. Harvell, P.F. Sale, A.J. Edwards, K. Caldeira, N. Knowlton, C.M. Eakin, R. Iglesias-Prieto, N. Muthiga, R.H. Bradbury, A. Dubi, M.E. Hatziolos, Coral reefs under rapid climate change and ocean acidification, Science 318 (2007) 1737–1742. http://dx.doi.org/10.1126/science.1152509. [8] P.W. Sammarco, K.B. Strychar, Responses to high seawater temperatures in zooxanthellae Octocorals, PLoS ONE 8 (2) (2013) e54989. http://dx.doi.org/10.1371/journal.pone. 0054989. [9] T.P. Hughes, A.H. Baird, D.R. Bellwood, M. Card, S.R. Connolly, C. Folke, R. Grosberg, O. Hoegh-Guldberg, J.B.C. Jackson, J. Kleypas, J.M. Lough, P. Marshall, M. Nystrm, S.R. Palumbi, J.M. Pandolfi, B. Rosen, J. Roughgarden, Climate change, human impacts, and the resilience of coral reefs, Science 301 (2003) 929–933. http://dx.doi.org/10.1126/science.1085046. [10] M.J.H. Van Oppen, J.C. Mieog, C.A. Sanchez, K.E. Fabricius, Diversity of algal endosymbionts (zooxanthellae) in octocorals: the roles of geography and host relationships, Molec. Ecol. 14 (2005) 2403–2417. http://dx.doi.org/10.1111/j.1365-294X.2005.02545.x. [11] P.W. Sammarco, K.B. Strychar, Effects of climate change on coral reefs: adaptation/exaptation in corals, evolution in zooxanthellae, and biogeographic shifts, Envtl. Bioindicators 4 (2009) 9–45. [12] E.P. Odum, Ecologia, Editora Guanabara, Rio de Janeiro, 1988. [13] R.H. Whittaker, Communities and Ecosystems, second ed., Macmillan, New York, 1975. [14] J.L. Harper, The Population Biology of Plants, Academic Press, New York, 1977. [15] P.L. Antonelli, R.H. Bradbury, Volterra–Hamilton Models in the Ecology and Evolution of Colonial Organisms, Series in Mathematical Biology and Medicine, World Scientific, Singapore, 1996. [16] G.E. Hutchinson, A note of the theory of competition between social species, Ecology 28 (1947) 319–321. [17] P.L. Antonelli, R.S. Ingarden, M. Matsumoto, The Theory of Sprays and Finsler Spaces with Applications in Physics and Biology, in: FTPH, 58, Springer/Kluwer, Dordrecht, 1993. [18] P.L. Antonelli, S.F. Rutz, C.E. Hirakawa, The mathematical theory of endosymbiosis I, Nonlinear Anal. RWA 12 (2011) 3238–3251. [19] P.L. Antonelli, S.F. Rutz, K.T. Fonseca, The mathematical theory of symbiogenesis II: models of the fungal fusion hypothesis, Nonlinear Anal. RWA (2012) http://dx.doi.org/10.1016/j.nonrwa.2012.01.005. [20] P.L. Antonelli, L. Bevilacqua, S.F. Rutz, Theories and models in symbiogenesis, Nonlinear Anal. RWA 4 (2003) 743–753. [21] G.D. Stanley Jr, P.K. Swart, Evolution of the coral-zooxanthenllae symbiosis during the Triassic: a geochemical approach—Paleobiology, 1995. [22] G.D. Stanley Jr, B. Van De Schootbrugge, The evolution of the coralalgal symbiosis, in: Coral Bleaching, Springer, Berlin Heidelberg, 2009, pp. 7–19.
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
73
[23] S.F. Rutz, Differential Geometric Modelling in Trophodynamics, in: Rubem Mondaini (Ed.), Proceedings of the First Brazilian Symposium on Mathematical and Computational Biology, 320 pp, ISBN 85-87922-26-2, E-papers, Rio de Janeiro, 2001. [24] P.L. Antonelli, P.W. Sammarco, Evolution of multiple complementary (secondary) metabolites, their synergism and their stability in colonial organisms, Open Sys. Inf. Dyn. 6 (1) (1999). [25] P.L. Antonelli, T.J. Zastawniak, Fundamentals of Finslerian Diffusion with Applications, in: FTPH, 101, Springer/Kluwer, Dordrecht, 1999. [26] P.L. Antonelli (Ed.), Handbook of Finsler Geometry, Vol. 2, Springer/Kluwer, Dordrecht, 2003. [27] P.L. Antonelli, B.C. Lackey (Eds.), The Theory of Finslerian Laplacians, in: Mathematics and its Applications, vol. 459, Springer/Kluwer series, Dordrecht, 1998. [28] M. Nagasawa, Segregation of a population in an environment, J. Math. Biol. 9 (2013) 213–235. [29] E. Nelson, Quantum Fluctuations, Princeton U. Press, 1985. [30] S.F. Rutz, R. Portugal, FINSLER: a computer algebra package for finsler geometry, Nonlinear Anal. 47 (2001) 6121–6134. and in [26] 2nd volume. [31] www.maplesoft.com. [32] J. Needham, Chemical Embriology, Vol. 3, Cambridge U. Press, 1931. [33] J. Huxley, Problems of Relative Growth, Dover Press, 1972. [34] D. Kosambi, Systems of differential equations of second order, Quart. J. Math. (Oxford Ser.) 6 (1935) 1–12. [35] E. Cartan, Observations sur le memoir precedent, Math. Zeitschrift 37 (1933) 619–622. [36] S.S. Chern, Sur la geometrie d’un systeme d’equations differentialles du second ordre, Bull. Sci. Math. 63 (1939) 206–212. [37] P.L. Antonelli, T.J. Zastawniak, Noise induced transitions in stochastic Volterra–Hamilton, open systems, Open Syst. Inf. Dyn. 4 (1997) 89–100.
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
A coral bleaching model Peter L. Antonelli a,∗,1 , Solange F. Rutz b , Paul W. Sammarco c , Kevin B. Strychar d a
Department of Mathematical Sciences, University of Alberta, Edmonton, AB, Canada
b
Federal University of Pernambuco - UFPE, Mathematics Department, Recife, PE, Brazil
c
Louisiana Universities Marine Consortium - LUMCON, Chauvin, LA, USA d Annis Water Resources Institute, Grand Valley State University, Muskegon, MI, USA
article
info
Article history: Received 25 May 2013 Accepted 14 September 2013
abstract Volterra–Hamilton systems theory is used to model coral bleaching. The algal and coral partners in the obligate symbiont organism exchange compounds, each one producing compounds benefiting the other. It is found that the production equations, based on a homogeneous cost functional, have unstable solutions which become stable in a constant temperature environment. Temperature spikes and global warming occur in the noisy background. The variety of clade distributions for the algal symbiont is shown to be a consequence of Nelson’s stochastic mechanics. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Coral bleaching is a result of the expulsion of symbiont algae, known as zooxanthellae, by their host organisms commonly called stony corals, or hermatypic scleractinians, the main reef builders [1]. Bleaching is a major problem globally because it threatens the existence of tropical reefs [2–4]. This process is called bleaching because of the loss of pigmentation contributed to the holobiont by the zooxanthellae, resulting in the white calcium carbonate exoskeleton of the coral host showing through the still intact but transparent coral tissue. Zooxanthellae exist in a large number of clades or variations of this alga and are found in the interstitial spaces of the coral tissue in all hermatypic or reef-building corals [5,6]. They are also found in soft corals, or octocorallians [7]. Zooxanthellar clades vary between coral species. There is always one clade of zooxanthellae in a given coral species in a given location and time and sometimes more than one clade. Corals can only survive and grow within a relatively narrow range of seawater temperatures and temperature increases predicted to occur in the world’s oceans could conceivably drive today’s reefs towards extinction [2]. This is the ecological problem addressed in this paper. Recent research has shown that the scleractinian host corals have adapted, over evolutionary time, to moderate fluctuations in seawater temperature, exhibiting a relatively broad range of temperature tolerances, while their symbiotic zooxanthellar partners have not [8]. Indeed, the zooxanthellae seem to be catching up, adapting more rapidly to temperature-based environmental fluctuations, which are increasing with time because of global warming [9–11]. The generation time (i.e., length of time it takes to reach age of first reproduction) for corals is on the order of years, while generation time for the zooxanthellae is on the order of days to weeks [11]. This enables the zooxanthellae to reproduce at a considerably greater rate than their coral hosts by possessing a much higher frequency of reproduction. A shorter generation
∗
Corresponding author. Tel.: +1 780 438 7886. E-mail addresses: [email protected] (P.L. Antonelli), [email protected] (S.F. Rutz), [email protected] (P.W. Sammarco), [email protected] (K.B. Strychar). 1 Visiting Professor at UFPE, Recife, PE, Brazil. 1468-1218/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.nonrwa.2013.09.006
66
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
time implies a better ability of these types of populations to adapt, merely on the basis of frequency of reproduction, all other things being equal. But, when generation time is compared to mutation rate in these organisms, it was determined that, on a per generation basis, zooxanthellae actually have a mutation rate lower than expected for an organism of their size [11]. This is counterintuitive, since mutations are one of the factors that drive adaptation and evolution. If this is the case, the question is raised, which of the two organisms comprising the holobiont, is actually responsible for most of the adapting? Could it be that the zooxanthellar mutation rate is sufficient to produce, on an evolutionary time scale, the many types of zooxanthellar clades now known to exist? We modelled these traits in order to understand and possibly predict survivorship and adaptation in corals under environmental conditions of rising seawater temperatures, known to be occurring at present. Our approach is both ecological and physiological in perspective. 2. Basics of our approach There are tens of algal clades which make up zooxanthellae, but for a coral species at a given reef location and time, only a few algal types at most are housed, and these can differ genetically and often in pigmentation. It has been suggested that the relatively short generation time of the zooxanthellae, together with natural selection, would provide a great number of variable zooxanthellar types, on an evolutionary time scale. Under normal environmental conditions, there would be a well-defined average value for their spatial density in the coral tissue, with a relatively small variance over time. We denote this density as, N 1 (t ), with t in a fixed interval [t0 , t1 ], in a single coral colony at a fixed location. Although there may be more than one type of zooxanthella living in symbiosis within a single coral, for simplicity in model building, we will assume here that there is just one. There is considerable variety and turnover of zooxanthellae in a coral—based upon repeated colonization. We presume this colonization process is controlled by interactions between the coral cells and the symbionts via an immune response, which is in turn genetically determined but may be influenced by the local, proximal, ecological and environmental conditions. This interaction can be described as in classical mathematical ecology using approximately constant coefficients in first order nonlinear differential equations [12–15]. There are two populations, one of coral cells, the other of algal cells, and these are taken to be weakly coupled which is expressed by coefficients being constant, or nearly so. This is in contrast to the strongly coupled case in which the coefficients actually depend explicitly on the x1 , x2 (amounts produced by, and exchanged with, each partner in a holobiont organism) and even on the ratio of the N [16,17]. Strong couplings do occur in nature and involve robustly stable, genetically programmed, chemical exchanges [18–20]. The term ‘‘robustly stable’’, as used here, means that the production stability measure, K , is positive and stays positive both before and after small perturbations of the x-values are performed in the expression for K . However, the bleaching situation addressed here seems not to be robustly stable, so there must be some small perturbations of x-values which render K non-positive. A specific zooxanthellar clade in residence in a specific coral tissue is not considered a permanent resident, despite the obligate relationship between the host and the symbiont. It might be referred to as a ‘‘fair weather symbiont’’ and a clade within a coral may change over time. This is far different from the obligate symbionts in, say, root fungi in vascular plants, or organelles, as mitochondria and chloroplasts in eukaryotes [18,19]. The coral acquire zooxanthellae as symbionts and these microalgae enable the corals to produce aragonitic calcium carbonate, the basic substrate of a coral reef. This ability is an evolutionary innovation which came about during the late Triassic period [21,22]. Before this, corals had no skeletons and did not form reefs (e.g. rugose corals). Nowadays, of course, each population in the holobiont gains from the association and some in ways naturally described as allometric. For example, the coral host #1, obtains photosynthates (i.e. sugars and oxygen) from the algae in an amount denoted x1 = b1 ln(m1 ) + c1 , while the algae, #2, receive nitrogenous wastes (e.g., CO2 and water) from the coral in the amount, x2 = b2 ln(m2 )+ c2 . Here, the b and the c are allometric constants and we say that the x are surrogates of biomass. Our model applies to all surrogates of biomass, not just this particular interpretation involving sugars and nitrogen. Further, we consider this obligate relationship as relatively weak and not linked or locked. Zooxanthellar turnover does certainly occur in their natural habitat. Indeed, corals expel zoox, mostly already dead—those which were not able to survive the rise in sea water temperature. One outstanding question is, does the holobiont organism replace expelled algae with a new viable zooxanthellar clade in time to permit continued reef calcium carbonate production? We believe that it does because it follows from the mathematical model below. We will assume holobiont production makes efficient use of available energy (ATP) for the system as a whole; so we will adjoin Volterra production equations, dxi /dt = k(i) N i to the ecological equation system of first order non-linear differential equations, which are to be chosen below [15,17]. This is the basic idea of analytical trophodynamics theory, which focuses on stability properties of the production processes [23,17]. Again, population #1 is the zooxanthellae and #2 is the coral species polyp population. (The parentheses on the lower index indicate that the equation is not summed from 1 to 2. However, the Einstein summation convention, summing on repeated upper and lower indices, is used throughout.) The parameters ki are positive proportionality constants and can be set equal to unity for convenience in this particular model. Since the combined system of equations is of second order, we can use the calculus of variations and implement simple regularity conditions which ensure physiological production processes
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
67
minimize production cost (ATP). We choose a cost function F which is positive and allow for explicit dependence on the produced substances, x1 , x2 , their production rates dx1 /dt, dx2 /dt, and on clock time, t. Moreover, the optimal growth of the two partner populations would, under ideal conditions, be in synchrony – strongly correlated – in order to ensure the spatial distribution of zooxanthellae within the coral cell population remains approximately constant. In reality, this may not be the case because of external factors which can alter the zooxanthellar populations (e.g., disease, ocean acidification, etc.). Here, we will assume that both N 1 and N 2 , are increasing initially at a fixed mean rate, λ. This rate is not to be confused with generation time. Moreover, the measure of joint production will be the cost itself. Thus, the infinitesimal relation ds = F (x1 , x2 , dx1 , dx2 , t , dt ) implies that small positive increments dx1 and dx2 , occurring in time interval, [t , t + dt ], will result in small positive increments ds, in joint production. We will further suppose that twice (or any positive multiple, C of) an increment, would multiply by a factor of two (or by a factor, C ), the joint production in that interval. This is known as the homogeneity assumption, and is biologically sound and has been used previously [24]. It should be noted that the parameter s is akin to size. It is independent of units used to describe it. Indeed, s is to be thought of as a geometric quantity. Thus, a smooth curve γ in the production space has a finite length s given by the integral of ds along γ from one point to another. The positive quantity s is the measure of total joint production along γ , for the entire holobiont population under study. From calculus of variations, the path integral of ds will depend on the smooth path chosen. Generally, we will consider those with dx1 /dt , dx2 /dt both positive and differentiate the integral of F with respect to the path and set it equal to zero. This preliminary will ensure efficiency of the physiologically based production process, because the usual technical assumption of regularity results in equations for the path of minimum cost, between (xio , to ) and (xi1 , t1 ), relative to the given cost functional F . These are the so-called Euler–Lagrange equations associated with F . However, as in the above discussion, these equations must have constant coefficients or nearly so, and this requirement strongly constrains the possible choices of F . There are just three possibilities according to the classification theorem of Antonelli and Matsumoto (see appendix, [25]). However, two of these cost functions are not appropriate, because, in one case, it leads to asymmetrical and biased interactions, like parasitic or commensal ones, and in the other, there are no interactions at all, only two independent logistic equations, one for each population. Only one of the three cost functions leads to an a priori unbiased treatment description of the two cell populations. Namely, ds = [(dx1 )2 + (dx2 )2 ]1/2 exp{(L2 + 1)αi xi + λt /2 + L arctan(dx1 /dx2 ) + ψ(x)}, with L as defined below and ψ(x) = (1/2)[ν1 .(x1 )2 + ν2 (x2 )2 ] [25]. Here, all coefficients denoted by the α , the ν and λ are positive constants. The corresponding Euler–Lagrange equations split naturally into Volterra’s production equations above, with ki = 1, ψ(x) = 0 and ecological equations dN 1 /dt = λN 1 − 2(α2 )N 1 N 2 − (α1 )[(N 1 )2 − (N 2 )2 ], dN 2 /dt = λN 2 − 2(α1 )N 1 N 2 − (α2 )[(N 2 )2 − (N 1 )2 ]. Of course, substituting the Volterra production equations back into these, results in the geodesic equations of the Finsler geometry automatically determined by the cost F . However, the path parameter here, is clock time t, not parameter s (see [15], p. 47). The α , along with λ, determine the conceptually deeper (geometric) structures of the cost F . Note that if N 2 is set to zero in equation #1, then a logistic differential equation is obtained for N 1 , so the algal population will decrease to a (relatively) small size (we assume small carrying capacity K1 ) and will die out from noisy environmental fluctuations alone. Likewise, a logistic is obtained from equation #2, if N 1 = 0 and N 2 will decrease to (relatively) small size (we assume a small carrying capacity K2 ) and will be wiped out by noisy fluctuations. Note that we take, α1 = λ/K1 and α2 = λ/K2 , where K1 and K2 are the (relatively) small carrying capacities for independently living populations of algal and coral cells, respectively. This smallness expresses their relative inability to survive on their own. But, they are partners in a symbiosis that is obligate and this ensures their survival. Remark. Do not confuse ki with Ki or K . They are very different! There is a unique steady-state in the positive (N 1 , N 2 ) quadrant when ψ(x) = 0. This is easy to compute and we leave it for the reader. The constant L = (α2 − α1 )/(α1 + α2 ) may be zero, negative or positive, but the production process is always unstable, since ψ(x) = 0 implies, K = 0, from the K -formula given below. Therefore, it follows there are paths in production space, close to start with, but which diverge eventually. Positive curvature, on the other hand, ensures that all initially close paths remain close. Negative or zero curvature forces instability so that at least some holobiont populations have divergent production curves. If K is strictly less than zero throughout production space, all paths must eventually diverge. Divergent paths, however, are unreliable from an evolutionary viewpoint since reproduction would be unstable, it being just a special case of production in our model, since the x-variables are surrogates of biomass with an assumed allometric relationship between reproductive and total biomass. Remark. It is of interest to write the parameters L and L2 + 1 in terms of the carrying capacities K1 and K2 . They are, L = (K1 − K2 )/(K1 + K2 ) and L2 + 1 = 2(K12 + K22 )/(K1 + K2 )2 . It is obvious that if the zooxanthellae carrying capacity K1 is less than K2 of the coral, then L is negative in the expression for K , while the reverse is true when K2 < K1 . Some standard references for Finsler geometry are [26,17,15,27,25].
68
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
3. Motivation for the model For the unstable case, we have ψ(x) = 0. Yet, certain negative constant terms representing temperature, added to the right-hand sides of the ecological equations, convert the production process to one that is stable. We prove this below using KCC-theory [20]. But, before we do this, we clarify several points: Firstly, if we had positive curvature K (and therefore production stability), then the holobiont would be highly integrated, physiologically. The reason is that positivity of curvature K can only be achieved by including ψ(x) = (1/2).[ν1 (x1 )2 + ν2 (x2 )2 ] in the exponential function and, as a consequence, the coefficients in the ecological equations must explicitly depend on xi , and the interactions, consequently, would exhibit strong coupling. Although there is solid evidence that a chemical exchange between the two components of this symbiotic system is part of an evolutionary innovation for the holobiont, the chemical exchanges occurring are weak compared to the highly integrated robustly stable exchanges mentioned above. Nevertheless, they are obligate, because the hermatypic corals cannot survive without their zooxanthellae, and conversely. Secondly, because there are several clades of zooxanthellae available for the replacement of one currently held by a given species of coral, when and if a seawater temperature spike occurs, and since the (joint) production in the holobiont is considered weak, any stability must be due of the influence of the ambient environment. This is demonstrated below. Third and finally, the model constructed here predicts that (joint) production stability in the holobiont can be achieved, but only over evolutionary time-scales and resulting from physiological changes and integration of the resulting chemical processes into the holobiont genome. Furthermore, the mathematical model predicts that the environmentally forced turnover of zooxanthellar types within a single intact coral colony, does take place, and is caused by rising seawater temperatures and temperature spikes. This is a conclusion from the mathematical model resulting from environmental noise having been added to the model according to the stochastic theory for conservative potential driven systems [28,29]. 4. Constant temperature case The stability measure K is actually the Gauss–Berwald curvature K of the ‘‘production space’’ for the holobiont. The stability scalar K has been computed, via S.F. Rutz’s software package FINSLER [30] based on MAPLE [31]. She expresses it as K = −2(ν1 + ν2 )(L2 + 1) exp{[2φ(x) + 2L arctan(N 1 /N 2 )]},
φ(x) = ψ(x) + 2(L2 + 1)[αi xi ], L = (α2 − α1 )/(α1 + α2 ), L2 + 1 = [(α1 )2 + (α2 )2 ]/(α1 + α2 )2 . It is significant that L and K together are sufficient to uniquely specify the Finsler geometry of the production surface [26,17]. Let us discuss the N-ratio effect on the curvature K . If ν1 and ν2 were positive and, if in addition, N 1 /N 2 was small, it follows that arctan(N 1 /N 2 ) is positive and close to zero which then implies that K , being negative, would be closer to zero. That is, the process of production, although unstable, progresses towards stability. Thus, the reapportionment of the total energy available to the holobiont to include relative population size, indicated by N 1 /N 2 , would help to stabilize production. The importance of this for holobiont evolution is clear since stable production is a necessary condition for stable reproduction. Let us explain further. It is standard procedure in Volterra–Hamilton systems theory (i.e., analytical trophodynamics) to assume that x1 , x2 are surrogates of total accumulated biomass [23,17]. This allows allometric quantification of the products, for example, seed production in plants, fertilized egg production in animals, production of chemicals during growth or even morphological characters. Thus, reproductive materials would be in allometric proportion to the total biomass produced. This is familiar to biologists in the guise of linear regression in log–log plots of phenotypic characters, going back to J. Needham and J. Huxley in the 1930s [32,33]. We argued above that incorporation of negative quadratic terms in φ(x), which necessarily resulted in stability of production, is not well founded. On similar grounds, inclusion of positive versions of φ(x) is not well founded, either. Therefore, we must take both ν1 and ν2 to be negligibly small or zero. This brings us right back to the system with the α restricted to be more or less constant, and for which small perturbations in x-values may alter the stability. 5. KCC-theory proof of environmentally induced stability Let (x1 , . . . , xn ) = (x), (dx1 /dt , . . . , dxn /dt ) = (dx/dt ) = (˙x), and t be 2n + 1 coordinates in an open connected subset Ω of the Euclidean (2n + 1)-dimensional space Rn × Rn × R1 . Suppose that we have d2 xi dt 2
+ g i (x, x˙ , t ) = 0,
i = 1, . . . , n
(1)
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
69
for which each g i is C oo in a neighbourhood of initial conditions ((x)0 , (˙x)0 , t0 ) ∈ Ω . The intrinsic geometric properties of (1) under non-singular transformations of the type
x¯ i = f i (x1 , . . . , xn ), t¯ = t ,
i = 1 , . . . , n,
(2)
are given by the five KCC-differential invariants, named after D. Kosambi [34], E. Cartan [35], and S.S. Chern [36], given below. Let us first define the KCC-covariant differential of a contravariant vector field ξ i (x) on Ω by
Dξ i
=
dt
dξ i dt
1
+ g;ir ξ r ,
(3)
2
where the semi-colon indicates partial differentiation with respect to x˙ r , and we made use of the Einstein summation convention on repeated indices. Using (3), Eq. (1) becomes
D x˙ i dt
= ϵi =
1 2
g;ir x˙ r − g i ,
(4)
defining the first KCC-invariant of (1), the contravariant vector field on Ω , ϵ i , which represents an ‘external force’. Varying trajectories xi (t ) of (1) into nearby ones according to x¯ i (t ) = xi (t ) + ξ i (t )η,
(5)
where η denotes a parameter, with |η| small and ξ (t ) the components of some contravariant vector defined along x = x (t ), we get, substituting (5) into (1) and taking the limit as η → 0, i
d2 ξ i dt 2
+ g;ir
dξ r dt
i
+ g,ir ξ r = 0,
i
(6)
where the comma indicates partial differentiation with respect to xr . Using the KCC-covariant differentiation (3) we can re-express this as
D 2ξ i dt 2
= Pri ξ r ,
(7)
where
Pji = −g,ij −
1 2
g r g;ir ;j +
1 2
x˙ r g,ir ;j +
1 4
g;ir g;rj +
1 ∂ 2 ∂t
g;ij .
(8)
The tensor Pji is the second KCC-invariant of (1). The third, fourth and fifth invariants are:
1 i i i Rjk = (Pj;k − Pk;j ) , 3 i Bjkl = Rjki ;l , i Djkl = g;ij;k;l .
(9)
The main result of KCC-theory is the following: Two systems of the form (1) on Ω are equivalent relative to (2) if and only if the five KCC-invariants are equivalent. In particular, there exists coordinates (¯x) for which the g i (¯x, x˙¯ , t ) all vanish if and only if all KCC-invariants are zero. The tensor D vanishes if and only if g i is quadratic in x˙ , in the case when the first KCC-invariant vanishes. Let us now introduce the notion of a n-dimensional Finsler space as a manifold where, given a coordinate system (x) and a curve xi = xi (t ), the norm of a tangent vector x˙ i to the curve at each point P on xi (t ) is given by the positive metric function F , |˙xi | = F (x, x˙ ), where F is positively homogeneous of degree 1 in x˙ i . From F , a metric tensor is defined as gij (x, x˙ ) = (∂ 2 F 2 /∂ x˙ i ∂ x˙ j )/2, which must be regular in an open region of the tangent bundle, the collection of all tangent vectors to the manifold, and which excludes the origins. The use of the calculus of variations for F leads to (1) with g i (x, x˙ , t ) = γjki (x, x˙ ) x˙ j x˙ k , where the γjki are the Levi-Cività symbols for the Finsler metric tensor gij (x, x˙ ). Berwald’s Gaussian curvature K for two-dimensional Finsler spaces is defined from his famous formula [17] i Rjk = F K mi (lj mk − lk mj ),
(10)
i where Rjk is given by the first equation in (9), li = x˙ i /F is the unit vector in the x˙ i direction, and mi the unique (up
to orientation) unit vector perpendicular to li . Lowering the index on mi via the metric tensor gives mi , which satisfies F (x, m) = gij (x, x˙ ) mi mj ≡ mi mi = 1. If our curvature K is bigger than zero everywhere, then trajectories oscillate back and forth, crossing the reference trajectory. In this case, we say (1) is Jacobi stable. If K ≤ 0 everywhere, trajectories diverge and system (1) is Jacobi unstable [26,15]. This notion of stability is a Lyapunov notion, but it is a whole trajectory concept.
70
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
We consider the RHS of the production equations above with L not necessarily zero. The RHS will be of the form,
−∇ i (V (x)), where V depends only on x-values and the gradient is of the form: g ij (x, dx/dt )∂j with gij = 1/2(∂˙i )(∂˙j )F 2 . The resulting equations are Euler–Lagrange equations for the Lagrangian
L = (1/2)gij (dxi /dt )(dxj /dt ) − V (x). The so-called metric tensor gij depends on the ratios N 1 /N 2 = dx1 /dx2 and is globally defined on the symmetric tensor product T (M ) ⊗ T (M ) with T (M ) denoting the tangent bundle of positive conical regions of the tangent planes, all with origin removed. This is why the non-zero rate condition must be implemented (see [26] vol. 1, p. 570). The deviation curvature tensor for the conservative dynamics of L is found to be
−Pji = Rjki .yk + ∇ i (Vj (x)), V (x) = a/2 exp[2αi xi ],
Vj (x) = ∂j V (x),
g ki = exp[−2αm xm ]h(y)ki
∇ k (Vj (x)) = a.h(y)ki [2(αi )(αj ) − (αr )Grij ] := a.hki Qij . Here, hki depends only on yj = dxj /ds and is positive definite on the tangent cones, and furthermore, Giii = αi ; Giij = αj , if i ̸= j; Gijj = −αi , if i ̸= j. These three-index quantities are symmetric in the two lower indices. It is seen that for all four choices of i, j Qij = [(α1 )2 + (α2 )2 ], so that the gradient of V (x) is positive definite. Of course, Rijk is vanishing as follows from the formula relating it to K . We can now see that the deviation vector ξ i , which is orthogonal to dxi /ds of γ , is oscillating back and forth across the reference trajectory γ in production space and therefore the production process in terms of the allometric x-variables is stable. 6. Nelson’s stochastic mechanics and zooxanthellar clades Define on the Euclidean 2-plane, R2 , the potential V = (a/2) exp[αi xi ]. In fact,
−a(α i ) = −g ij (x)(∂j )(V (x)) which are the right-hand sides of the Euler–Lagrange equations for Lagrangian L = (1/2)gij yi yj − V (x), where yi = dxi /ds and g ij = exp{−2αk xk }. δ ij is the inverse of gij , a Riemannian metric tensor. The Euler–Lagrange equations are dx1 /ds = y1 ,
dy1 /ds = −2(α2 )y1 y2 + (α1 )[(y2 )2 − (y1 )2 ] − a(α1 ),
dx2 /ds = y2 ,
dy2 /ds = −2(α1 )y1 y2 + (α2 )[(y1 )2 − (y2 )2 ] − a(α2 ).
The stochastic theory used here was developed by Edward Nelson of Princeton University and his students for a stochastic approach to derivation of Schrödinger’s equation in classical quantum mechanics. The paths are nowhere differentiable but are continuous and Markov. This theory is mathematically difficult, but has been much studied by mathematical physicists [29]. However, a surprising biological application to social behaviour was presented in [28]. Here, we show how our bleaching model can be transformed into the double harmonic oscillator considered in [28]. A discrete spectrum eigenvalue problem, with known solutions, is obtained. Transforming coordinates, spatial distributions for successively higher total energy are obtained with each a sort of ‘‘mountain range’’ with the lowest energy level yielding a single Gaussian distribution in the phenotypic character space spanned by our surrogates of biomass variables (see appendix, [15]. There are nodal lines defined by setting coordinates equal to zero, and the enclosed regions are isolated from each other in the sense that crossing a nodal boundary, by any sample path starting inside, has probability zero. However, the processes are ergodic within these compartments [15,29,28]. Theorem. There is a non-singular coordinate transformation from x1 , x2 to polar coordinates r , θ given by r (x1 , x2 ) = exp(α1 )x1 + (α2 )x2 ,
θ (x1 , x2 ) = (α2 )x1 − (α1 )x2 , so that our environmentally perturbed system becomes d2 r (s)/ds2 − r [dθ (s)/ds]2 = −kr (s), d2 θ (s)/ds2 + (2/r (s))[dr (s)/ds][dθ (s)/ds] = 0, where k = a[(a1 )2 + (a2 )2 ] and s(t ) = A − B exp{−(λ)t } are convenient definitions of constants and path parameter s in terms of clock time t. This is the polar form of the double harmonic oscillator.
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
71
Fig. 1. Density ρ0 (x, y).
The Newton–Nelson formula [26, page 269] is equivalent to Schrödinger’s equations for the two dimensional harmonic oscillator, but Planck’s constant is replaced by ν 2 , and the stationary solution ψ satisfies
−
ν2 2
∆+V
ψ = λψ,
√ ρ . According to well-known sources, the discrete eigenvalues satisfy √ λn = ν 2 a (n + 1), n = 1, 2, . . . .
where ψ =
The stationary solution for the smallest eigenvalue λ0 is
λ0 = C0 exp −
1 2
a
ν2
r
2
.
The corresponding density function is a single ‘‘bump’’. For the second eigenvalue λ1 there are two independent stationary solutions (see Figs. 1 and 2 and [37] for more graphs), and for the third, λ2 , there are three, and so forth. Finally, the densities, as pictured, are given by squaring the stationary solutions. These densities describe chemical morphology of the holobiont. The upshot of this method of noise modelling is that the total energy (ATP plus increments produced by the gradient of an environmental potential V ) yields distributions (the diagrams above) for the holobiont population, depending on the total energy level, with the larger energies yielding more peaks (more ‘‘mountains in the mountain range’’). Most of the holobiont variety, as discussed above, is understood to be the various clades of the symbiont algae. Furthermore, when a random (seawater temperature) spike event occurs in the environment of a holobiont state situated near a compartmental boundary, the occupied state can continuously move towards the nearby boundary and pass through to a new distribution. Yet, all compartments at a given total energy level are isolated from each other in Nelson stochastic theory. By allowing higher and higher total energy levels, all holobiont types occurring in nature can, in principle, be represented. A random seawater temperature spike can force communication between the isolated compartments. As time goes on, more and more ‘‘mountains’’ are occupied providing the necessary rich assortment of types. To be precise, the parameter a is assumed to increase or decrease by jumps of positive or negative magnitude δ , at random times, the number J+ (S ) of jumps by +δ and the number J− (S ) of jumps by −δ up to time s, being Poisson processes independent of each other and of the background noise in the environment. The requisite computer simulations were published previously [37]. The standard approach to adding noise to Finslerian production equations has been developed for several distinct two-dimensional geodesic systems [25]. However, the Nelson theory for conservative systems has not been developed and not for lack of trying! Fortunately, we are able in this present model to use the Riemannian version successfully. Detailed proofs have appeared in several places, in particular, [15, appendix], [25,26], or [37].
72
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
Fig. 2. Density ρ01 (x, y).
Acknowledgement FACEPE grant BFP-0043-1.01/11. References [1] C.E. Birkeland, The Life and Death of Coral Reefs, Springer, New York, 1997, p. 536. [2] O. Hoegh-Guldberg, Climate change, coral bleaching, and the future of the world’s coral reefs, Mar. Freshw. Res. 50 (1999) 839–866. [3] M.J.H. Van Oppen, J.C. Mieog, C.A. Sanchez, K.E. Fabricius, Diversity of algal endosymbionts (zooxanthellae) in octocorals: the roles of geography and host relationships, Molec. Ecol. 14 (2005) 2403–2417. http://dx.doi.org/10.1111/j.1365-294X.2005.02545.x. [4] M.J.H. Van Oppen, J.M. Lough (Eds.), Coral Bleaching: Patterns, Processes, Causes, and Consequences, in: Ecol. Studies, vol. 205, Springer, New York, 2009, p. 178. [5] R. Rowan, Review—diversity and ecology of zooxanthellae on coral reefs, J. Phycol. 34 (1998) 407–417. http://dx.doi.org/10.1046/j.15298817.1998.340407.x. [6] A.C. Baker, Flexibility and specificity in coral-algal symbiosis: diversity, ecology, and biogeography of symbiodinium, Ann. Rev. Ecol. Evol. Systematics 34 (2003) 661–689. http://www.jstor.org/stable/30033790. [7] O. Hoegh-Guldberg, P.J. Mumby, A.J. Hooten, R.S. Steneck, P. Greenfield, E. Gomez, C.D. Harvell, P.F. Sale, A.J. Edwards, K. Caldeira, N. Knowlton, C.M. Eakin, R. Iglesias-Prieto, N. Muthiga, R.H. Bradbury, A. Dubi, M.E. Hatziolos, Coral reefs under rapid climate change and ocean acidification, Science 318 (2007) 1737–1742. http://dx.doi.org/10.1126/science.1152509. [8] P.W. Sammarco, K.B. Strychar, Responses to high seawater temperatures in zooxanthellae Octocorals, PLoS ONE 8 (2) (2013) e54989. http://dx.doi.org/10.1371/journal.pone. 0054989. [9] T.P. Hughes, A.H. Baird, D.R. Bellwood, M. Card, S.R. Connolly, C. Folke, R. Grosberg, O. Hoegh-Guldberg, J.B.C. Jackson, J. Kleypas, J.M. Lough, P. Marshall, M. Nystrm, S.R. Palumbi, J.M. Pandolfi, B. Rosen, J. Roughgarden, Climate change, human impacts, and the resilience of coral reefs, Science 301 (2003) 929–933. http://dx.doi.org/10.1126/science.1085046. [10] M.J.H. Van Oppen, J.C. Mieog, C.A. Sanchez, K.E. Fabricius, Diversity of algal endosymbionts (zooxanthellae) in octocorals: the roles of geography and host relationships, Molec. Ecol. 14 (2005) 2403–2417. http://dx.doi.org/10.1111/j.1365-294X.2005.02545.x. [11] P.W. Sammarco, K.B. Strychar, Effects of climate change on coral reefs: adaptation/exaptation in corals, evolution in zooxanthellae, and biogeographic shifts, Envtl. Bioindicators 4 (2009) 9–45. [12] E.P. Odum, Ecologia, Editora Guanabara, Rio de Janeiro, 1988. [13] R.H. Whittaker, Communities and Ecosystems, second ed., Macmillan, New York, 1975. [14] J.L. Harper, The Population Biology of Plants, Academic Press, New York, 1977. [15] P.L. Antonelli, R.H. Bradbury, Volterra–Hamilton Models in the Ecology and Evolution of Colonial Organisms, Series in Mathematical Biology and Medicine, World Scientific, Singapore, 1996. [16] G.E. Hutchinson, A note of the theory of competition between social species, Ecology 28 (1947) 319–321. [17] P.L. Antonelli, R.S. Ingarden, M. Matsumoto, The Theory of Sprays and Finsler Spaces with Applications in Physics and Biology, in: FTPH, 58, Springer/Kluwer, Dordrecht, 1993. [18] P.L. Antonelli, S.F. Rutz, C.E. Hirakawa, The mathematical theory of endosymbiosis I, Nonlinear Anal. RWA 12 (2011) 3238–3251. [19] P.L. Antonelli, S.F. Rutz, K.T. Fonseca, The mathematical theory of symbiogenesis II: models of the fungal fusion hypothesis, Nonlinear Anal. RWA (2012) http://dx.doi.org/10.1016/j.nonrwa.2012.01.005. [20] P.L. Antonelli, L. Bevilacqua, S.F. Rutz, Theories and models in symbiogenesis, Nonlinear Anal. RWA 4 (2003) 743–753. [21] G.D. Stanley Jr, P.K. Swart, Evolution of the coral-zooxanthenllae symbiosis during the Triassic: a geochemical approach—Paleobiology, 1995. [22] G.D. Stanley Jr, B. Van De Schootbrugge, The evolution of the coralalgal symbiosis, in: Coral Bleaching, Springer, Berlin Heidelberg, 2009, pp. 7–19.
P.L. Antonelli et al. / Nonlinear Analysis: Real World Applications 16 (2014) 65–73
73
[23] S.F. Rutz, Differential Geometric Modelling in Trophodynamics, in: Rubem Mondaini (Ed.), Proceedings of the First Brazilian Symposium on Mathematical and Computational Biology, 320 pp, ISBN 85-87922-26-2, E-papers, Rio de Janeiro, 2001. [24] P.L. Antonelli, P.W. Sammarco, Evolution of multiple complementary (secondary) metabolites, their synergism and their stability in colonial organisms, Open Sys. Inf. Dyn. 6 (1) (1999). [25] P.L. Antonelli, T.J. Zastawniak, Fundamentals of Finslerian Diffusion with Applications, in: FTPH, 101, Springer/Kluwer, Dordrecht, 1999. [26] P.L. Antonelli (Ed.), Handbook of Finsler Geometry, Vol. 2, Springer/Kluwer, Dordrecht, 2003. [27] P.L. Antonelli, B.C. Lackey (Eds.), The Theory of Finslerian Laplacians, in: Mathematics and its Applications, vol. 459, Springer/Kluwer series, Dordrecht, 1998. [28] M. Nagasawa, Segregation of a population in an environment, J. Math. Biol. 9 (2013) 213–235. [29] E. Nelson, Quantum Fluctuations, Princeton U. Press, 1985. [30] S.F. Rutz, R. Portugal, FINSLER: a computer algebra package for finsler geometry, Nonlinear Anal. 47 (2001) 6121–6134. and in [26] 2nd volume. [31] www.maplesoft.com. [32] J. Needham, Chemical Embriology, Vol. 3, Cambridge U. Press, 1931. [33] J. Huxley, Problems of Relative Growth, Dover Press, 1972. [34] D. Kosambi, Systems of differential equations of second order, Quart. J. Math. (Oxford Ser.) 6 (1935) 1–12. [35] E. Cartan, Observations sur le memoir precedent, Math. Zeitschrift 37 (1933) 619–622. [36] S.S. Chern, Sur la geometrie d’un systeme d’equations differentialles du second ordre, Bull. Sci. Math. 63 (1939) 206–212. [37] P.L. Antonelli, T.J. Zastawniak, Noise induced transitions in stochastic Volterra–Hamilton, open systems, Open Syst. Inf. Dyn. 4 (1997) 89–100.