Interactive comment on “Differences and implications in biogeochemistry from maximizing entropy production locally versus globally” by

In this manuscript, the author propose a simply two reaction, coupled model that represents autotrophy and metabolism. Then using Monte Carlo approach to investigate the model domain space and demonstrate some important results. 1) that there are many micropath solutions that produce the same macropath state and 2) global optimum does not correspond to local optimums. Overall, this paper is an elegant and important work that contributes to our understanding of the thermodynamic development of simple ecological systems. It could be published with minor revisions.


Introduction
There is a long history of research that attempts to understand and model ecosystems from a goal-based perspective dating back to at least Lotka (1922), who argued that ecosystems organize to maximize power.Many of the proposed theories (e.g., Morowitz, 1968;Schr ödinger, 1944;Odum and Pinkerton, 1955;Margalef, 1968;Prigogine and Nicolis, 1971;Schneider and Kay, 1994;Bejan, 2007;Jorgensen et al., Figures Back Close Full are inspired by or derived from thermodynamics, because ecosystems are comprised by numerous agents, and there is the desire to understand their collective action, which is an inherently thermodynamic-like approach.Thermodynamics appears particularly relevant for understanding biogeochemistry as it is largely under microbial control, or more aptly molecular machines (Falkowski et al., 2008), that span the realm between pure chemistry and biology.While significant effort and progress has been made in understanding organismal metabolism in an ecosystem context, it is rather surprising that we still lack an agreed upon theory that can explain and quantitatively predict ecosystem biogeochemistry that is independent of the constituent organisms.Perhaps a single theory does not exist, and ecosystem biogeochemistry depends completely or substantially on the nuances of the constituent organisms that comprise an ecosystem.However, at the appropriate scale and perspective it appears likely that systems organize in a predictable manner independent of species composition, as the planet's biogeochemical processes have remained relative stable over hundreds of millions of years but the organisms that comprise ecosystems have changed significantly over this period.Consider primary producers in terrestrial versus marine ecosystems.While there is a definitive morphological discontinuity between marine systems that are composed of planktonic algae and terrestrial ecosystems comprised of trees, shrubs or grasses, the distinction and discontinuity is removed when considering solar energy acquisition and carbon dioxide fixation, which both do similarly.Because of the large degrees of freedom, it is likely impossible to predict the nature of the organisms that evolve to facilitate energy and mass transformations in one ecosystem versus the next or how their populations vary over time, but at the level of free energy dissipation and biogeochemistry, the problem may be more tractable.The problem is similar to predicting weather versus climate, where the latter can be predicted over long time scales at the expense of details necessary for accurate weather prediction whose accuracy decays rapidly with time.Our objective is to understand biogeochemistry at the level of generic biological structure catalyzing chemical reactions.Understanding energy and Figures

Back Close
Full mass flow catalyzed by living systems is particularly relevant for understanding how the biosphere will respond to global change, and predicting biogeochemical responses is essential before any geoengineering projects could be responsibly implemented.Many of the thermodynamically inspired ecosystem theories are more similar than different (Fath et al., 2001;Jorgensen, 1994), but the theory of maximum entropy production (MEP) has made great progress in understanding how biotic and abiotic systems may organize under nonequilibrium steady state conditions when sufficient degrees of freedom exist.It appears Paltridge (1975) was the first to apply MEP as an objective function for modeling global heat transport, but it was not until the theoretical support provided by Dewar (2003) that research in MEP theory and applications garnered greater interest, particularly in Earth systems science (Lorenz et al., 2001;Kleidon et al., 2003;Kleidon and Lorenz, 2005a,b;Dyke and Kleidon, 2010).While Dewar's initial work has received some criticisms, there have been several other approaches that arrive at the same nonequilibrium steady state result of MEP (see Niven, 2009 and references therein).The MEP principle states that systems will organize to maximize the rate of entropy production, which is an appealing extension to the second law of thermodynamics for nonequilibrium systems.In essence, nonequilibrium systems will attempt to approach equilibrium via the fastest possible pathway, which can include the organization of complex structures if they facilitate entropy production.However, none of the current MEP theories incorporate information content of the system, such as that contained within an organism's genome, even though several of the MEP analyses are derived from Jaynes' (2003) maximum entropy (MaxEnt) formulation that is information based.Furthermore, current MEP theory suggests that maximizing entropy production locally will maximize entropy production over the system domain, at least for flow-controlled systems (Niven, 2009) In this manuscript we will assume that biological systems organize to maximize entropy production, which is synonymous with maximizing the dissipation rate of Gibbs free energy for systems under constant temperature and pressure.Living organisms are viewed here as mere catalysts facilitating autocatalytic reactions for the dissipation of free energy.We will also continue our main hypothesis (Vallino, 2010), that the acquisition of Shannon information (Shannon, 1948), or more precisely useful information (Adami, 2002;Adami et al., 2000), via evolution facilitates the production of entropy when averaged over time.It is the ability of living systems to store information that allows them to out-compete abiotic systems in entropy production under appropriate conditions.Here, we will examine how information may also facilitate entropy production averaged over space.For our model system, we will consider a simple two-box model where entropy production in each box is governed by two reactions that determined the rate of biological structure (i.e., catalyst) synthesis and degradation.

Model system
To examine how total entropy production depends on either local or global entropy production optimization, we use a simple two-box model reminiscent of ocean biogeochemistry models developed to capture processes in the photic and aphotic zones of oceans (Fig. 1) (Ianson and Allen, 2002).Model focus is placed on dissipation of free energy by synthesis and destruction of biological structures ( 9Entropy of mixing is readily calculated between boxes and boundaries based on flux and 11 chemical potential obtained from concentration differences between boxes (Meysman and12 Bruers 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit 13 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17 coefficient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or 18 boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21

) 22
Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 23 efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6)(7)(8)(9)(10)(11)(12). 25 1 and 9 Entropy of mixing is readily calculated between boxes and boundaries based on flux and 11 chemical potential obtained from concentration differences between boxes (Meysman and12 Bruers 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit 13 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17 coefficient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or 18 boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21

) 22
Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 23 efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6-12).25 2 ) from and to chemical constituents (CH 2 O, NH 3 , O 2 and H 2 CO 3 ) within the two layers.The photic zone is typically energy replete but resource limited, while the aphotic zone is its complement: energy limited but resource replete.Our model captures these two differing states; however, we replace energy input via photosynthetic active radiation with energy input via the diffusion of chemical potential from a fixed upper boundary condition into the surface layer box.Nitrogen in the form of ammonia is allowed to diffuse into the All constituents are allowed to freely diffuse between boxes [1] and [2], but only CH 2 O, O 2 and H 2 CO 3 are permitted across the upper boundary, while only NH 3 can diffuse across the lower boundary (Fig. 1).
In our simplified model system, we consider only two autocatalytic chemical reactions occur within each box, as constrained by the following stoichiometric equations: , divided by temperature gives the entropy production per unit area for mixing between box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing can be calculated from Eqs. (6-12).
, divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21 Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 23 efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6)(7)(8)(9)(10)(11)(12). 25 (1) The product of flux and chemical potential differences, s the entropy production per unit area for mixing between , as follows (J m -2 d -1 ºK -1 ) 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth , entropy production associated with reactions and mixing .
erwise it is 0. The product of flux and chemical potential differences, rature gives the entropy production per unit area for mixing between l species k, as follows (J m -2 d -1 ºK -1 ) 1 and S 2 along with the growth each box, entropy production associated with reactions and mixing qs.(6-12).
2 ) oundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, , divided by temperature gives the entropy production per unit area for mixing between   (6)(7)(8)(9)(10)(11)(12). 25 (2) In reaction r 1 , biological structure, ere  is a characteristic length scale for diffusion (2.5×10 m), D is the diffusion [j]; otherwise it is 0. The product of flux and chemical potential differences, j i k , , divided by temperature gives the entropy production per unit area for mixing between 1 and S 2 along with the growth iciencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing n be calculated from Eqs. (6-12).
1 , catalyzes its own synthesis from sugar (unit carbon basis) and available ammonia thereby fulfilling the primary producer role.Reaction r 2 represents heterotrophic organisms that consume biological structure synthe-10 sized from Eq. (1), 9 ) therwise it is 0. The product of flux and chemical potential differences, perature gives the entropy production per unit area for mixing between ical species k, as follows (J m -2 d -1 ºK -1 ) 1 , but also functions in a cannibalistic mode by consuming 9 ( )   (6)(7)(8)(9)(10)(11)(12). 25 2 as well.While cannibalism could be replaced by inclusion of a sufficient number of higher trophic levels for closure, this would lead to greater degrees of model freedom that are not important for this study, but are more typical closure techniques observed in real environments (Edwards and Yool, 2000).To minimize parameterization, we assume 15 that biological structure is degraded indiscriminately by 9 Bruers 2007, Kondepudi and Prigogine 1998).For our simple model, we only permit 13 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 ( ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17 boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6-12).25 2 , so that δ i in Eq. ( 2) are given by, 9 lculated between boxes and boundaries based on flux and m concentration differences between boxes (Meysman and rigogine 1998).For our simple model, we only permit xes and boundaries, so a material flux of chemical species k x or boundary [j] is given by (mmol m -2 d -1 ), ) gth scale for diffusion (2.5×10 -5 m), D is the diffusion it is 0. The product of flux and chemical potential differences, gives the entropy production per unit area for mixing between ies k, as follows (J m -2 d -1 ºK -1 ) , NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth box, entropy production associated with reactions and mixing -12). ) cale for diffusion (2.5×10 -5 m), D is the diffusion The product of flux and chemical potential differences, s the entropy production per unit area for mixing between , as follows (J m -2 d -1 ºK -1 ) ) wise it is 0. The product of flux and chemical potential differences, ture gives the entropy production per unit area for mixing between species k, as follows (J m -2 d -1 ºK -1 ) ) [j]; otherwise it is 0. The product of flux and chemical potential differences, y temperature gives the entropy production per unit area for mixing between chemical species k, as follows (J m -2 d -1 ºK -1 ) ) [j]; otherwise it is 0. The product of flux and chemical potential differences, y temperature gives the entropy production per unit area for mixing between chemical species k, as follows (J m -2 d -1 ºK -1 ) where C tivities (Alberty 2003).
lated between boxes and boundaries based on flux and concentration differences between boxes (Meysman and ogine 1998).For our simple model, we only permit and boundaries, so a material flux of chemical species k r boundary [j] is given by (mmol m -2 d -1 ), ) scale for diffusion (2.5×10 -5 m), D is the diffusion compound k is allowed between box or s 0. The product of flux and chemical potential differences, es the entropy production per unit area for mixing between k, as follows (J m -2 d -1 ºK -1 ) H 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth , entropy production associated with reactions and mixing ).
i is the concentration (mmol m −3 ) of biological structure i .In both reactions, ρ is the N:C ratio of biological structure, terms of ( 9) and (10) instead of activities (Alberty 2003).
Entropy of mixing is readily calculated between boxes and boundaries based on flux and chemical potential obtained from concentration differences between boxes (Meysman andBruers 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit diffusive transport between boxes and boundaries, so a material flux of chemical species k from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion  , and ε i is the growth efficiency for biolog-20 ical structure synthesis.Note, as growth efficiency approaches zero, both reactions represent the complete oxidation of reduced carbon to CO 2 and H 2 O.Typically, Holling-type kinetics (Holling, 1965) are used to describe reaction rates for Eqs. ( 1) and (2) based on a limiting substrate or substrates.However, kinetic parameters (e.g., µ Max , K S ) in Holling-type expressions depend on community composition, Figures

Back Close
Full so Holling kinetic equations need to be reparameterized as community composition changes.Since an underlying hypothesis of the MEP principle when applied to living systems is that they will evolve and organize to achieve a MEP state (Vallino, 2010), we cannot a prior set parameters in a kinetic model as we do not know the nature of the community composition nor the corresponding kinetic parameter values that will lead to a MEP state.Kinetic parameters need to be dynamic to reflect community composition changes.To achieve this objective, we have developed a novel kinetic expression that can capture reaction kinetics in a manner consistent with community compositional changes.The form of our kinetic expression takes the hyperbolic form of the Monod equation (Monod, 1949), where substrate specific up-take rate is given by (d −1 ), and specific growth rate readily follows as (d −1 ), In Eq. ( 4), C j is the concentration (mmol m −3 ) of one or more growth limiting substrates, such as NH 3 , ε i is the growth efficiency as used in Eqs. ( 1) and (2), and ν * and κ * are universal parameters that remain fixed regardless of community composition.Equation (4) differs from the Monod equation in several important ways.The effective half saturation "constant" (κ * ε 4 i ) depends on growth efficiency.The rational for this dependency is that nutrients at low concentration require free energy expenditure to transport them into the cell against a concentration gradient; consequently, organisms that have evolved to grow under low nutrients conditions (i.e., K-selection; Pianka, 1970) will, by thermodynamic necessity, grow at lower efficiencies.The (1 − ε 2 i ) term in Eq. (4) approximately accounts for thermodynamic constraints on kinetics, because net reaction rate must approach zero as reaction free energy tends toward zero (Jin and Bethke, 2003;Boudart, 1976) term in Eq. ( 4) is empirically motivated to account for reduced substrate uptake rates at low growth efficiency.The universal parameters have been chosen to qualitatively match observations of bacterial growth.Under oligotrophic conditions often found in ocean gyres', bacterial specific growth rate is typically 1-2 d −1 , with growth efficiencies of 10-20%, and substrate concentrations in the µM or lower range (Del Giorgio and Cole, 1998;Carlson et al., 1999).At the other extreme, bacteria grown under ideal laboratory conditions show specific growth rates as fast as 50 d −1 , with growth efficiencies around 50-60%, provided substrate concentrations are in the mM range (Bailey, 1977;Lendenmann and Egli, 1998).Equations ( 4) and ( 5) capture these approximate trophic extremes 10 with values of κ * of 5000 mmol m −3 and ν * of 350 d −1 (Fig. 2).The usefulness of the kinetic expressions (Eqs.4 and 5) is their ability to generate organismal growth kinetics expected under oligotrophic conditions to extreme eutrophic conditions by varying only growth efficiency, ε (Fig. 2).Based on Eq. ( 4), the rate of Eqs. ( 1) and ( 2) are expressed by (mmol m −3 d −1 ), 15 Bruers 2007, Kondepudi and Prigogine 1998).For our simple model, we only permit 13 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 () where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17 [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21

) 22
Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 23 efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6-12).25 f mixing is readily calculated between boxes and boundaries based on flux and potential obtained from concentration differences between boxes (Meysman and 07, Kondepudi and Prigogine 1998).For our simple model, we only permit transport between boxes and boundaries, so a material flux of chemical species k or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) s a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, ided by temperature gives the entropy production per unit area for mixing between ing is readily calculated between boxes and boundaries based on flux and tial obtained from concentration differences between boxes (Meysman and ondepudi and Prigogine 1998).For our simple model, we only permit port between boxes and boundaries, so a material flux of chemical species k ) haracteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion ; otherwise it is 0. The product of flux and chemical potential differences, by temperature gives the entropy production per unit area for mixing between r chemical species k, as follows (J m -2 d -1 ºK -1 ) 1 and S 2 along with the growth and 2 ε in each box, entropy production associated with reactions and mixing ed from Eqs. (6-12).
Entropy of mixing is readily calculated between boxes and boundaries based on flux and chemical potential obtained from concentration differences between boxes (Meysman andBruers 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit diffusive transport between boxes and boundaries, so a material flux of chemical species k from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion [j]; otherwise it is 0. The product of flux and chemical potential differences, , divided by temperature gives the entropy production per unit area for mixing between box

Free energy and entropy production
For our two box model, entropy production occurs via destruction of chemical potential via Eqs.( 1) and ( 2), and by entropy of mixing associated with relaxing chemical gradients between boxes [1] and [2] and their associated boundaries.For chemical reactions, entropy production per unit area in a given box with depth h (m) is readily determined under constant pressure and temperature from the reaction Gibbs free energy change, ∆G r , reaction rate, r, and absolute temperature, T , as given by Eu (1992, 131-141 While Gibbs free energy of reaction is usually straight forward to calculated for a given reaction stoichiometry, both Eqs. ( 1) and ( 2) require some discussion because both involve synthesis or decomposition of living biomass, or more generally, what we refer to as biological structure.
Living organisms are often thought to be very low entropy structures because they are highly order systems.However, this common misconception results from the incorrect association of thermodynamic entropy with order or disorder in a system.In general, system order has little to do with thermodynamic entropy (Kozliak and Lambert, 2005;Lambert 1999), as the latter is purely related to dispersion of energy to lower energy states, which may or may not be less ordered.Interestingly, both theoretical calculations and empirical data show that bacteria and yeast have higher absolute entropies (at 298.15 • K) than glucose (Battley, 1999a(Battley, ,b, 2003)).Biological structure has a higher entropy of formation than "simple" growth substrates, and it can be shown that the standard Gibbs free energy of reaction for synthesizing bacteria or yeast from glucose, ammonia, phosphate and sulfate, without CO 2 production, is slightly less than zero, so is a reaction that can occur spontaneously (Vallino, 2010).Hence, it is thermodynamically possible for Eq. ( 1) to proceed without additional free energy input with ε 1 = 1.However, as mentioned above, reaction rates are also constrained by thermodynamics

ESDD Introduction Conclusions References
Tables Figures

Back Close
Full as Gibbs free energy of reaction approaches zero.For a reaction to proceed at high rate it must proceed irreversibly; consequently, organisms whose growth efficiency is less than 100% can grow faster than those operating at very high efficiencies.The diminishing returns of high growth efficiency is embed in our growth kinetics equation, as plotting Eq. ( 5) as a function of ε exhibits an optimum, so that maximum specific The Gibbs free energy of reaction for Eqs. ( 1) and ( 2) can be expressed as a linear combination of biosynthesis and CH 2 O oxidation, as given by (J mmol −1 ), [j]; otherwise it is 0. The product of flux and chemical potential differences, , , divided by temperature gives the entropy production per unit area for mixing between en concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth iencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing be calculated from Eqs. (6-12).
where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17   ( ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion [j]; otherwise it is 0. The product of flux and chemical potential differences, , divided by temperature gives the entropy production per unit area for mixing between box Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing can be calculated from Eqs. (6)(7)(8)(9)(10)(11)(12). ) )  ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6-12).25 where R is the ideal gas constant (in J mmol • Ox is the standard Gibbs free energy of CH 2 O oxidation by O 2 to H 2 CO 3 (−498.4J mmol −1 ) and ∆G • 9 chemical potential obtained from concentration differences between boxes (Meysman and12 Bruers 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit 13 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 ( )  is the standard Gibbs free energy for biological structure synthesis from CH 2 O and NH 3 , which is set to zero for our study because it is negligibly small (Vallino, 2010).We use the approach of Alberty (2003Alberty ( , 2006) ) to calculate the standard Gibbs free energy 15 of reaction, ∆G • Ox , which accounts for proton dissociation equilibria between chemical species (CO 2 + H 2 O ↔ H 2 CO 3 ↔ H + + HCO − 3 , etc.) at a pH of 8.1 and temperature of 293.15 • K. Ionic strength (I S = 0.7 M) is used to estimate activity coefficients, so that concentrations can be used in the logarithmic correction terms of Eqs. ( 9) and ( 10) instead of activities (Alberty, 2003).(Meysman and Bruers, 2007;Kondepudi and Prigogine, 1998).For our simple model, we only permit diffusive transport between boxes and boundaries, so a material flux of chemical species k from box or boundary [i ] to box or boundary [j ] is given by (mmol m −2 d −1 ), where is a characteristic length scale for diffusion (2.5 Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , 9 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 ( ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17 coefficient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or 18 boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21

Transport equations, optimization and numerical routines
To obtain changes in chemical constituents over time for a given set of boundary conditions, a mass balance model is constructed based on in-and-out fluxes and production rates by reactions r 1 and r 2 for each constituent k in boxes [1] and [2], which takes the general form, (1) and ( 2), and r is a 2 × 1 vector of reaction rates, [r 1 r 2 ] T , given by Eqs. ( 6) and ( 7).The full model equations are detailed in Appendix A for k = CH 2 O, O 2 , H 2 CO 3 , NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth ntropy production associated with reactions and mixing 1 , NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth ox, entropy production associated with reactions and mixing 2).
2 in the two boxes, and parameter values used for simulations is given in Table 1.

5
We have intentionally designed our biogeochemistry model to contain only a few transport related parameters (D, , h) and effectively only two adjustable biological parameters per box, giving a total of four adjustable parameters for the two box model: , and ε 2 [2] that maximize entropy production can be determined by numerical analysis, and via Eqs.( 6) and ( 7) the biogeochemistry associated with the MEP state is defined.For a zero dimensional system, a variation of this approach has been used to determine how biogeochemistry develops over time 15 (Vallino, 2010); however, for a system that has one or more spatial dimensions, a choice needs to be made regarding local versus global optimization of entropy production.For the two-box model, entropy production can be asynchronously maximized locally (i.e., in each box), as given by, or globally (i.e, across both boxes synchronously) as, where both optimizations are bounded within the hypercube, 0 ≤ ε i [j ] ≤ 1.Of course, local maximization, as given by Eq. ( 14), requires an iterative (i.e., asynchronous) approach, because changing the conditions in one box alters the optimal solution for Figures

Back Close
Full the other.In contrast, global optimization given by Eq. ( 15) is mathematically well posed, but the parameter space is larger (4-D vs. 2-D in this case).Neither Eq. ( 14) nor Eq. ( 15) include the entropy of mixing contribution given by Eq. ( 12), as it is uncertain how it should be partitioned, but as will be shown, it is a small contribution compared to entropy production from the reaction terms, so neglecting it in the optimization is acceptable.
While optimization of Eqs. ( 14) and ( 15) can be conducted for the transient problem (Eq.13) over a specified time interval (Vallino, 2010), we are primarily interested in how spatial optimization alters MEP solutions.Consequently, we are only interested in steady state (SS) solutions to Eq. ( 13).Typically, Newton's method, or a variation thereof, is used to solve for SS solutions to Eq. ( 13); however, employing several variations to Newton's method including a line search method with trust region (Bain, 1993) was found to fail in some subspaces of 0 ≤ ε i [j ] ≤ 1.An approach based on interval arithmetic (Kearfott and Novoa, 1990;Kearfott, 1996) was also developed, but proved to be too computationally burdensome.The most robust approach found was to integrate the ODE equations (Eq.13) forward in time using a high precision block implicit method (Brugnano and Magherini, 2004) from a specified initial condition ( Numerical solution to Eq. ( 15) was obtained by using VTDIRECT (He et al., 2009), which employs a parallel version of a Lipschitzian direct search algorithm (Jones et al., 1993) to find a function's global optimum without using function derivatives.We also used VTDIRECT to solve Eq. ( 14), but we implemented an iterative approach, where first box 3 Results and discussion

Model characteristics
For any specified point in ε i [j ]-space, a solution to Eq. ( 13) can be found.For instance, Fig. 3 illustrates the transient dynamics to steady state for the point ε 2 [2]) = (0.1, 0.05, 0.1, 0.05) with the parameters and initial and boundary conditions 5 used throughout this study (Tables 1 and 2, respectively).For this particular solution, SS entropy production dominates in box ), respectively.Due to the significant amount of free energy re-10 lease in oxidation of reduced organic compounds, entropy of mixing is a small fraction of the entropy production associated with Eqs. ( 1) and (2).To demonstrate this more generally, we randomly chose 100 000 points within the 4-D ε i [j ]-space and calculated the steady state solution to Eq. ( 13) and associated entropy terms, which shows that entropy of mixing is at most 7.6% of the entropy of reaction, and for most cases much 15 less than that (Fig. 4a).Only 310 points out of 100 310 random simulations did not attain a SS solution within 10 6 days.
An interesting aspect of the model is that the amount of nitrogen extracted from the sediments, or boundary [3], depends on community composition, because changing values of ε i [j ] dramatically affects the total standing SS nitrogen within boxes [1] 20 and [2], calculated as (mmol N m −2 ), ate activity coefficients, so that concentrations can be used in the logarithmic correction of ( 9) and ( 10) instead of activities (Alberty 2003).
py of mixing is readily calculated between boxes and boundaries based on flux and ical potential obtained from concentration differences between boxes (Meysman andrs 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit sive transport between boxes and boundaries, so a material flux of chemical species k box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) e  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion icient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or dary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, j , divided by temperature gives the entropy production per unit area for mixing between i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) n concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth encies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing e calculated from Eqs. (6-12).
estimate activity coefficients, so that concentrations can be used in the logarithmic correction 8 terms of ( 9) and ( 10) instead of activities (Alberty 2003).9 10 Entropy of mixing is readily calculated between boxes and boundaries based on flux and 11 chemical potential obtained from concentration differences between boxes (Meysman and12 Bruers 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit 13 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 () where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17 coefficient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or 18 boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21

) 22
Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 23 efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6-12).25 which is evident in the results from the 100 000 random ε i [j ]-point simulations (Fig. 4b).
For any given reaction entropy production defined by Eq. ( 8), numerous SS solutions can be found that differ in location within ε i [j ]-space and exhibit differences of up Introduction

Conclusions References
Tables Figures

Back Close
Full to four orders of magnitude in N Total (Fig. 4b).While the vast majority of solutions result in an N Total less than 10 4 mmol m −2 , 33 solutions produced high to extremely high N Total values.Examination of community kinetics that lead to high N Total values reveals an apparent necessary condition that the growth efficiency of Eq. ( 2), ε 2 [j ], in both boxes [1] and [2] must be less than 0.02 (Fig. 5).Inspection of the transient 5 solution reveals that high 9 d by temperature gives the entropy production per unit area for mixing between for chemical species k, as follows (J m -2 d -1 ºK -1 ) ntrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 1 ε and 2 ε in each box, entropy production associated with reactions and mixing ated from Eqs. (6-12).
1 concentration develops due to extremely slow growth rates of 9 s 0. The product of flux and chemical potential differences, es the entropy production per unit area for mixing between k, as follows (J m -2 d -1 ºK -1 )     ε in each box, entropy production associated with reactions and mixing be calculated from Eqs. (6-12).
2 .Inspection of total entropy production by Eqs. ( 1) and ( 2) in the random 100 000 10 ε i [j ]-point simulations reveals that the top 17 SS solutions have total reaction entropy production distributed over a narrow range from 620.11 to 621.15 J m −2 d −1 • K −1 ; however, the corresponding ε i [j ] points have a much broader distribution over ε i [j ]-space (Table 3), which produce significant differences in transient dynamics and SS standing stocks (Fig. 6).Most significant are differences in SS standing stocks for biological 15 structures 9 es and boundaries, so a material flux of chemical species k or boundary [j] is given by (mmol m -2 d -1 ), ) th scale for diffusion (2.5×10 -5 m), D is the diffusion it is 0. The product of flux and chemical potential differences, ives the entropy production per unit area for mixing between es k, as follows (J m -2 d -1 ºK -1 ) NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth ox, entropy production associated with reactions and mixing 12).
1 and 9 ween boxes and boundaries, so a material flux of chemical species k [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) istic length scale for diffusion (2.5×10 -5 m), D is the diffusion therwise it is 0. The product of flux and chemical potential differences, perature gives the entropy production per unit area for mixing between cal species k, as follows (J m -2 d -1 ºK -1 )

Global versus local MEP solutions
The global MEP solution obtained by maximizing reaction associated entropy production, Eq. ( 8), summed over boxes [1] and [2] by simultaneously varying ε as given by Eq. ( 15) is 621.4).Because of the small 4-D parameter space of our model, this solution was almost located in the random 100 000 point simulations ( . This result indicates that the global optimum does not correspond to local optimums.Entropy productions associated with mixing, Eq. ( 12), summed over all relevant state variables for the globally-optimized solution are 18.6, 0.06 and 0.0 J m −2 d −1 • K −1 for fluxes F 0,1 , F 1,2 and F 2,3 , respectively.While entropy production from mixing does contribute to total entropy production, it is clear that it is small relative to the entropy of reaction terms.
If reaction associated entropy production is maximized locally as given by Eq. ( 14), then the total entropy production summed over both boxes is only 429.4 J m −2 d −1 • K −1 , with 224.4 and 204.9 J m −2 d −1 • K −1 being produced in boxes [1] and [2], respectively (Table 4).Examination of entropy production in boxes [1] and [2] in the vicinity of the locally-optimized solutions (Fig. 8) illustrates that the iterative procedure used to solve Eq. ( 14 the solution is stable.The corresponding entropy of mixing terms for fluxes F 0,1 , F 1,2 and F 2,3 are 3.21, 2.60 and 0.0 J m Comparing steady state solutions from the global versus the local optimizations shows that CH 2 O concentration in boxes [1] and [2] is near zero in the globallyoptimized solution, but only partially depleted in the locally-optimize solution (Table 5).

5
The inability to effectively oxidize CH 2 O in box [1] in the locally-optimized solution results from an insufficient development of biological structure ( 9boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21

) 22
Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 23 efficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6-12).25 2 ), as it is evident in Fig. 4b that a lower boundary on N Total exists for a given entropy production rate.Biological reaction rates are ultimately limited by the quantity of catalyst available that in turn is limited by elemental resources.The spatial model configuration 10 is such that state variables in box [1] control energy acquisition across boundary [0], while box [2] controls resource input across boundary [3].Limitations in either flux will result in lower entropy production rates.In the local optimization given by Eq. ( 14), maximizing entropy production in box [2] ultimately results in less N transported into the system from boundary [3].In the globally optimize solution, net synthesis of from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 ()   (6)(7)(8)(9)(10)(11)(12). 25 1 [2] 15 is critical in N acquisition, because when reaction r 1 exceeds r 2 in box [2] there is a net transport of NH 3 into the system that accumulates in 9 diffusive transport between boxes and boundaries, so a material flux of chemical species k 14 from box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), 15 ()   (6)(7)(8)(9)(10)(11)(12). 25 1 [2].In the locally optimize solution, entropy production can be increased by a large loop flow between reactions r 1 and r 2 , but this limits N acquisition and entropy production in box [1].These simulations clearly demonstrate that optimizing local entropy production limits whole system 20 resource and energy acquisition, which results in lower entropy production compared to the globally optimized solution that more effectively utilizes available resources to dissipate free energy.
We also investigated globally and locally optimized solutions with different values of the transport parameters (namely D, , h[1] and h [2]), but in each case the general con-25 clusion was supported; entropy production is greater when maximized globally rather than locally.A second variation of the model was also investigated that included two additional state variables, )   9) and ( 10) instead of activities (Alberty 2003).
ntropy of mixing is readily calculated between boxes and boundaries based on flux and hemical potential obtained from concentration differences between boxes (Meysman andruers 2007, Kondepudi andPrigogine 1998).For our simple model, we only permit iffusive transport between boxes and boundaries, so a material flux of chemical species k rom box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) here  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion oefficient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or oundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, j i k , µ , divided by temperature gives the entropy production per unit area for mixing between ox [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) for this model was that and S 2 along with the growth d 2 ε in each box, entropy production associated with reactions and mixing from Eqs. (6-12).
i optimized in box [j ] would be poorly adapted to conditions in box [k] and would be outcompeted, but could serve as substrate in reaction r 2 .In this alternate model, diffusion across boxes occurred between active and inactive variants of , H 2 CO 3 , S 1 and S 2 along with the growth py production associated with reactions and mixing i , as given by 9 ical species k, as follows (J m d ºK ) en concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth ciencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing be calculated from Eqs. (6-12).
Next, we examine how our approach and conclusions may be extrapolated to a more ecologically relevant context.

Broader ecological context
The model we have developed is obviously chosen to illustrate various aspects in the 10 application of MEP in spatially explicit situations, and does not represent any natural system per say.The ideas we explore here, based on the development and results of our modeling exercise, are intended to examine the usefulness of the MEP principle for understanding biogeochemical processes.The two main questions the MEP community must address are (1) do biological systems organize towards a MEP state and

15
(2) if so, is this knowledge useful for understanding and modeling biogeochemistry?This manuscript does not address question (1), as we implicitly assume it to be true; instead, we focus on question (2).Ecosystem processes and the resulting biogeochemistry are often viewed from an organismal perspective, because organisms are the autonomous parts that comprise

Conclusions References
Tables Figures

Back Close
Full MEP state under a given set of constraints (Dewar, 2003).Consequently, ecosystem biogeochemistry is not constrained to one particular set of species, but can be attained by an infinite number of complementary species sets, with the sets being interchangeable.By complementary we mean the community must be comprised of organisms capable of energy acquisition and element recycling, so that at steady state only free 5 energy is dissipated.The Earth's biosphere obviously operates in this mode, as there is no net accumulation or loss of biomass over time, but the conversion of electromagnetic radiation into longer wavelengths produces entropy.We have attempted to embody the idea of species composition fluidity in kinetic expressions Eqs. ( 4) and ( 5), where the choice of ε i defines the kinetics of the community.To implement MEP, dif-10 ferent modeling approaches need to be developed that are independent of species composition, which Eq. ( 4) is but one example.
The MEP approach also highlights the importance of resource acquisition for the construction of biological structures needed to dissipate free energy.Different community parameterizations lead to different levels of N acquisition and biological struc-15 ture concentrations with a minimum 9 m box or boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) here  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion efficient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or undary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, j i k , µ , divided by temperature gives the entropy production per unit area for mixing between iven concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth ficiencies 1 ε and 2 ε in each box, entropy production associated with reactions and mixing n be calculated from Eqs. (6-12).
T concentration required to achieve a given rate of entropy production (Fig. 4b).This perspective differs from the current paradigm of Liebig's Law of the Minimum, which focuses solely on elemental limitations at the organismal level.Because organism stoichiometry changes to match resource availability over evolutionary time scales (Elser et al., 2000(Elser et al., , 2007)), Liebig's law becomes 20 a weak constraint.From the MEP perspective, resources will be extracted based on catalyst compositions that achieve the highest rate of free energy dissipation as constrained by organic chemistry.To quote Lineweaver and Egan ( 2008) "This represents a paradigm shift from "we eat food" to "food has produced us to eat it"." The 17 highest entropy producing solutions in the 100 000 random-point simulations 25 illustrate that very different "food web" configurations, as defined by ε i [j ], can nevertheless dissipate free energy at statistically similar rates (Table 3).Each of these MEP solutions can be considered as meta-stable or as alternate stable states, because a perturbation of sufficient magnitude, or alteration of initial conditions, can lead the Figures

Back Close
Full system to a different MEP basin of attraction where different ecosystem processes are exhibited for the same boundary conditions, such as occurs in shifts between macroalgae and coral reefs (Mumby, 2009), in microbial communities (Price and Morin, 2004) and in many other systems (Schr öder et al., 2005).Likewise, altering system drivers or boundary conditions will change the ecosystem configuration necessary to maintain a MEP state, which would likely lead to regime shifts (Brock and Carpenter, 2010).All of the solutions in Table 3 are stable states, as they were obtained via integration of the ODEs (Appendix A) to steady state.Unlike most analyses based on a discrete selection of food web components with fixed parameter values, our approach based on Eqs. ( 4) and ( 5) selects from a continuum of possible food web configurations that could arise over long periods of evolution, but does not produce a statistically unique solution (Table 3).While we have not investigated stability aspects of the MEP solutions in this manuscript, the MEP perspective may be useful in such analyses.
The optimization criteria implicit in MEP allows parameter values to be determined for a given set of environmental conditions and drivers instead of being fit to experimental data (e.g., Table 4).This improves model robustness as parameter values can be determined for conditions where experimental data are lacking, or where model resolution is coarse compared to the underlying physics, such as in modeling global heat transport (Paltridge, 1975;Kleidon et al., 2003).The optimization approach also replaces knowledge on information content in a system.As stated in the introduction, genomic information allows biological systems to out-compete abiotic systems under some situations when averaged over time and/or space.However, it is not just Shannon information, but useful information (Adami et al., 2000) that is relevant; e.g., genes that allow sulfate to serve as an electron acceptor are useless if sulfate is not present in the environment.This context dependency of genomic information complicates decoding an ecosystem's metagenome.While great strides are being made at reading and decoding environmental genomic data (DeLong, 2009;Gianoulis et al., 2009), we are still far from using that information to predict ecosystem biogeochemistry (Keller, 2005;Frazier et al., 2003).By assuming all known metabolic capabilities are Figures present in a system, MEP-based optimization can determine which pathways may be expressed and when (Vallino, 2010); furthermore, genomic information can be incorporated into the optimization as constraints when genomic surveys indicate the lack of certain metabolic capabilities.Hence, MEP theory can be used to link genomic content to free energy dissipation.
An interesting result from our MEP-based modeling exercise is that entropy production can be increased when optimized globally versus locally.However, to achieve high entropy production requires that the system organize such that certain locations within the model domain function "sub-optimally", as evident in Fig. 7b.That is, free energy dissipation in box [2] is lower than possible, but this allows greater entropy production in box [1], which more than offsets the lower entropy production in box [2].This "altruistic behavior" is inconsistent with standard Darwinian evolution, which places optimization at the level of the individual; however, our results are consistent with ideas on evolution of cooperation and multilevel selection theory (Bastolla et al., 2009;Clutton-Brock, 2009;Goodnight, 1990;Hillesland and Stahl, 2010;Nowak, 2006;Traulsen and Nowak, 2006;Wilson and Wilson, 2008).For instance, bacterial colonies are large spatial structures (∼cm) compared to the bacteria (∼µm) that comprise them, yet it is well established that bacteria produce quorum sensing compounds that cause bacteria on the periphery of a colony to express different genes than those in the colony's center (Camilli and Bassler, 2006;Keller and Surette, 2006;Shapiro, 1998).Communication between plants via emission of volatile organic compounds has also been demonstrated to improve plant defenses over spatial scales greater than the individual (Heil and Karban, 2010).Below ground, a single fungus can attain immense size and occupy up to 1000 ha with hyphae (Ferguson et al., 2003).Given that mycorrhizae can integrate with plant roots (Maherali and Klironomos, 2007;Vandenkoornhuyse et al., 2007;Whitfield, 2007) suggests some level of below ground communication occurs across entire forests.From a physics perspective, even an individual multicellular organism represents a large spatial structure that is 10 clear that information exchange over large spatial scales exists in nature, so it is conceivable that ecosystems could coordinate over space to maximize energy acquisition and dissipation.If ecosystem information exchange occurs over large spatial scales, and if this communication facilitates increased free energy dissipation, then the domain of a model needs to be chosen judiciously.Specifically, the model domain is defined by the spatial extent of information exchange, largely dictated by chemical signaling.In certain systems, such as microbial mats (van Gemerden, 1993), it seems reasonable to model these systems as spatially connected systems, while other systems, such as the surface and deep ocean, it is uncertain to what extent, if any, they are informationally connected over space.Furthermore, it is possible that entropy production from local optimization is statistically similar to that from global optimization, so there may be little gained from spatial communication, and local optimization (e.g., Follows et al., 2007) would be sufficient to describe system dynamics.We do not know the answer to this question, as all current biogeochemical models depend solely on local conditions and are Markovian in nature; more research is needed.

Conclusions
The maximum entropy production principle (Paltridge, 1975;Dewar, 2003) proposes that nonequilibrium abiotic or biotic systems with many degrees of freedom will organize towards a state of maximum entropy production, which is synonymous with maximizing free energy dissipation for chemical systems.The usefulness of MEP theory for understanding and modeling biogeochemistry is still under debate, and numerous applied and theoretical challenges need to be address before MEP becomes a common tool for solving problems in biogeochemistry.In previous work we have demonstrate how the MEP principle can be applied to biological systems under transient conditions provided entropy production is integrated over time (Vallino, 2010).In this manuscript we have examined how biogeochemical predictions differ if entropy Figures

Back Close
Full production is maximized locally versus globally for systems involving spatial dimensions.Biological systems explore biogeochemical reaction space primarily via changes in community composition and gene expression; consequently, we have developed a simple kinetic expression, given by Eq. ( 4), to capture changes in community composition that govern biogeochemistry by using a single parameter, ε i , that varies be-5 tween 0 and 1. Investigating MEP in a simple two-box biogeochemistry model involving two chemical reactions has revealed that globally optimized solutions generate higher entropy production rates than locally optimized solutions (Table 4), and there exists many alternate steady state solutions that produce statistically similar rates of entropy (Table 3).Results from globally optimized MEP solutions support hypotheses regard-Figures
entropy production per unit area for mixing between ollows (J m -2 d -1 ºK -1 ) 2 , H 2 CO 3 , S 1 and S 2 along with the growth opy production associated with reactions and mixing

2
[1] ]; otherwise it is 0. The product of flux and chemical potential differences, temperature gives the entropy production per unit area for mixing between emical species k, as follows (J m -2 d -1 ºK -1 ) ns of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 2 ε in each box, entropy production associated with reactions and mixing rom Eqs.(6-12).
For box [2] constituents, the equations are, between boxes and boundaries based on flux and ntration differences between boxes (Meysman and 1998).For our simple model, we only permit oundaries, so a material flux of chemical species k dary [j] is given by (mmol m -2 d -1 ), for diffusion (2.5×10 -5 m), D is the diffusion s 1 if flux of compound k is allowed between box or he product of flux and chemical potential differences, entropy production per unit area for mixing between ollows (J m -2 d -1 ºK -1 ) 2 , H 2 CO 3 , S 1 and S 2 along with the growth opy production associated with reactions and mixing ) ristic length scale for diffusion (2.5×10 -5 m), D is the diffusion otherwise it is 0. The product of flux and chemical potential differences, perature gives the entropy production per unit area for mixing between ical species k, as follows (J m -2 d -1 ºK -1 ) ncentrations can be used in the logarithmic correction s (Alberty 2003).
between boxes and boundaries based on flux and ntration differences between boxes (Meysman and 1998).For our simple model, we only permit oundaries, so a material flux of chemical species k dary [j] is given by (mmol m -2 d -1 ), for diffusion (2.5×10 -5 m), D is the diffusion s 1 if flux of compound k is allowed between box or he product of flux and chemical potential differences, entropy production per unit area for mixing between ollows (J m -2 d -1 ºK -1 ) 2 , H 2 CO 3 , S 1 and S 2 along with the growth opy production associated with reactions and mixing ficients, so that concentrations can be used in the logarithmic correction instead of activities (Alberty 2003).
readily calculated between boxes and boundaries based on flux and tained from concentration differences between boxes (Meysman and udi and Prigogine 1998).For our simple model, we only permit tween boxes and boundaries, so a material flux of chemical species k y [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) ristic length scale for diffusion (2.5×10 -5 m), D is the diffusion otherwise it is 0. The product of flux and chemical potential differences, perature gives the entropy production per unit area for mixing between ical species k, as follows (J m -2 d -1 ºK -1 ) 7, Kondepudi and Prigogine 1998).For our simple model, we only permit ransport between boxes and boundaries, so a material flux of chemical species k r boundary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( ) s a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion t (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, ided by temperature gives the entropy production per unit area for mixing between [j] for chemical species k, as follows (J m -2 d -1 ºK -1 )  (Alberty 2003).ted between boxes and boundaries based on flux and ncentration differences between boxes (Meysman and ine 1998).For our simple model, we only permit nd boundaries, so a material flux of chemical species k boundary [j] is given by (mmol m -2 d -1 ), ) ale for diffusion (2.5×10 -5 m), D is the diffusion ) k is 1 if flux of compound k is allowed between box or 0. The product of flux and chemical potential differences, the entropy production per unit area for mixing between as follows (J m -2 d -1 ºK -1 ) ) where  is a characteristic length scale for diffusion (2.5×10 -5 m), D is the diffusion 17 coefficient (1×10 -4 m 2 d -1 ) and ) ( , k j i β is 1 if flux of compound k is allowed between box or 18 boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19 , divided by temperature gives the entropy production per unit area for mixing between 20 box [i] to [j] for chemical species k, as follows (J m -2 d -1 ºK -1 ) 21 )  Full  3. Values in 4-D ε i [j ]-space that correspond to the highest entropy production (J m −2 d −1 • K −1 ) by reactions (r 1 ) and (r 2 ) summed over boxes [1] and [2] from the 100 000 random ε i [j ]-point simulations.Full  Full  9 obtained from concentration differences between boxes (Meysman and epudi and Prigogine 1998).For our simple model, we only permit between boxes and boundaries, so a material flux of chemical species k ary [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( )    quilibria between chemical species (CO 2 + H 2 O ↔ H 2 CO 3 ↔ H + + of 8.1 and temperature of 293.15ºK.Ionic strength (I S =0.7 M) is used to ficients, so that concentrations can be used in the logarithmic correction instead of activities (Alberty 2003).
readily calculated between boxes and boundaries based on flux and tained from concentration differences between boxes (Meysman and udi and Prigogine 1998).For our simple model, we only permit tween boxes and boundaries, so a material flux of chemical species k y [i] to box or boundary [j] is given by (mmol m -2 d -1 ), ( )  14).
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Because a MEP-based model can be implemented with either local or global optimization, this manuscript investigates if the type of spatial optimization affects the solution to a simple, MEP-based biogeochemistry modelDiscussion Paper | Discussion Paper | Discussion Paper | bottom box [2] from sediments, which fixes the lower boundary condition [3].Hence, energy flows into the system via box [1] and resources (nitrogen) flow in via box [2].Discussion Paper | Discussion Paper | Discussion Paper |

ε
in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs.
and boundaries based on flux and oncentration differences between boxes (Meysman and gine 1998).For our simple model, we only permit nd boundaries, so a material flux of chemical species k boundary [j] is given by (mmol m boxes and boundaries based on flux and d from concentration differences between boxes (Meysman and nd Prigogine 1998).For our simple model, we only permit n boxes and boundaries, so a material flux of chemical species k o box or boundary [j] is given by (mmol m

ε
in each box, entropy production associated with reactions and mixing from Eqs. (6-12).

ε
in each box, entropy production associated with reactions and mixing can be calculated from Eqs. (6-12).
Discussion Paper | Discussion Paper | Discussion Paper | , which is equivalent to a growth efficiency of 1.The leading ε 2 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

9
from box or boundary[i]  to box or boundary[j]  is given by (mmol m -2 d -1),

ε
in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs. (6-12).25

20
Entropy of mixing is readily calculated between boxes and boundaries based on flux and chemical potential obtained from concentration differences between boxes Figures Back Close Full Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

ε
in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs.
) where the brackets, [i ], following a variable denote the box-boundary location with [0] being the upper boundary and [3] being the lower boundary (Fig. 1).Elements of Discussion Paper | Discussion Paper | Discussion Paper | the 1 × 2 vector Λ k are the stoichiometric coefficients for compound k obtained from Eqs.
and ε 2[2].Because choosing a value for ε i [j ] effectively selects the kinetic characteristics of the community in box [j ], we can examine how entropy 10 production changes as a function of community composition.More importantly for this study, the values of ε 1 Discussion Paper | Discussion Paper | Discussion Paper | [1] was optimized for given values of ε 1 [2] and ε 2 [2], then box [2] was optimized given ε 1 [1] and ε 2 [1] from the previous optimization of box [1]; this iteration proceeded until ε[j ] − ε[j ] ≤10 −9 for each box, where is the Euclidian norm and the vector ε ( ε = [ ε1 ε2 ] T ) is the value of ε from the previous iteration.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) has indeed obtained a solution that simultaneously maximizes entropy production in both boxes, because d S r /d t in box [1] cannot be improved with the given values of ε 1 [2] and ε 2 [2], nor can d S r /d t in box [2] be improved given ε 1 [1] and ε 2 [1]; Discussion Paper | Discussion Paper | Discussion Paper |

ε
in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs.

ε
in each box, entropy production associated with reactions and mixing 24 can be calculated from Eqs.
flux of compound k is allowed between box or nd [j]; otherwise it is 0. The product of flux and chemical potential differences, d by temperature gives the entropy production per unit area for mixing between or chemical species k, as follows (J m -

ε
in each box, entropy production associated with reactions and mixing an be calculated from Eqs. (6-12).Discussion Paper | Discussion Paper | Discussion Paper | to[j]  for chemical species k, as follows (J m d ºK ) 21 of the model also produced similar results, where global optimiza-5 tion resulted in total entropy production of 620.9 J m −2 d −1 • K −1 , but local optimization generated entropy at a rate of only 415.0 J m Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3 to 10 6 times larger than the cells that comprise it, and surely those cells are spatially coordinated.Consequently, it is Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | boxes and boundaries based on flux and tained from concentration differences between boxes(Meysman and   udi and Prigogine 1998).For our simple model, we only permit tween boxes and boundaries, so a material flux of chemical species k y [i] to box or boundary[j]  is given by (mmol m -2 d -1 ), ) of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth 2 in each box, entropy production associated with reactions and mixing m Eqs.(6-12). 1 2 ) of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) 22Given concentrations of CH 2 O, NH 3 , O 2 , H 2 CO 3 , S 1 and S 2 along with the growth boundary[i]  to box or boundary[j]  is given by (mmol m d ), 15 Discussion Paper | Discussion Paper | Discussion Paper | Table Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .
Fig. 1.Schematic of the two box model with associated state variables and boundary conditions.All fluxes between boxes and boundaries are governed by diffustion (see Eq. 11) and only CH 2 O, O 2 and H 2 CO 3 are exchanged accross boundary [0] and box [1].Similarly, only NH 3 is allowed across box [2] and boundary [3].Boundary concentrations are held at fix values throughout the simulation (see Table2).

Fig. 3 .
Fig. 1.Schematic of the two box model with associated state variables and boundary conditions.All fluxes between boxes and boundaries are governed by diffustion (see Eq. 11) and only CH 2 O, O 2 and H 2 CO 3 are exchanged accross boundary [0] and box [1].Similarly, only NH 3 is allowed across box [2] and boundary [3].Boundary concentrations are held at fix values throughout the simulation (see Table2).

Fig. 5 .Fig. 6 .
Fig. 5. Points (triangles and squares) in 4-D ε i [j ]-space that resulted in total steady state N standing stock (N Total ) greater than 10 4 mmol m −2 .Triangles mark locations and magnitude of N standing stock for ε 1 and ε 2 in box [1], while squares mark locations and N standing stock for ε 1 and ε 2 in box [2].Also shown as small points are a subset of the ε 1 and ε 2 in box [1] (green) and box [2] (red) from the 100 000 randomly selected ε i [j ]-point simulations that have N Total < 10 4 mmol m −2 .Note, the y-axis range from 0 to 0.021 captures all points that meet the minimum N Total requirement of 10 4 mmol m −2 .
compound k is allowed between box or 18 boundary [i] and [j]; otherwise it is 0. The product of flux and chemical potential differences, 19

Table 2
T d C/d t ≤ ζ , where C is a vector of concentrations and ζ was set to 10 −8 mmol m −3 d −1 .If a steady state solution failed to be obtained after 10 6 days of integration, the point was flagged and removed from the search, but this rarely occurred.
, and entropies associated with mixing between boundary [0] and box [1], box [1] and box [2], and box [2] and boundary [3] are 14.6, 0.436 and 0.00 Table 3).To examine the solution in the vicinity of the global maximum, we plot d S r /d t in box [1] for any given value of ε 1 [1] and ε 2 [1] while holding ε 1 [2] and ε 2 [2] fixed at their global optimum values, and vice versa for box [2] (Fig. 7).It is clear from this plot that d S r /d t in box [1] of 619.6 J m −2 d −1 • K −1 is near the d S r /d t maximum in box [1] given ε 1 [2] and ε 2 [2] (Fig. 7a, light gray point); however, it is also evident that d S r /d t in box [2] that contributes 1.92 J m −2 d −1 • K −1 (Fig.7b, light gray point) to the global solution is far removed from the maximum d S r /d t possible in box [2] that could be attained given ε 1 [1] and ε 2 [1], which is approximately 236

Table 1 .
Model fixed parameter values.

Table 5 .
Steady state concentration (mmol m −3 ) of state variables at maximum entropy production associated with local and global solutions.