Photorespiration and Rate Synchronization in a Phototroph-Heterotroph Microbial Consortium

Abstract

The process of oxygenic photosynthesis is robust and ubiquitous, relying centrally on input of light, carbon dioxide, and water, which in many environments are all abundantly available, and from which are produced, principally, oxygen and reduced organic carbon. However, photosynthetic machinery can be conflicted by the simultaneous presence of carbon dioxide and oxygen through a process sometimes called photorespiration. We present here a model of phototrophy, including competition for RuBisCO binding sites between oxygen and carbon dioxide, in a chemostat-based microbial population. The model connects to the idea of metabolic pathways to track carbon and degree of reduction through the system. We find decomposition of kinetics into elementary flux modes a mathematically natural way to study synchronization of mismatched rates of photon input and chemostat turnover. In the single species case, though total biomass is reduced by photorespiration, protection from excess light exposures and its consequences (oxidative and redox stress) may result. We also find the possibility that a consortium of phototrophs with heterotrophs can recycle photorespiration byproduct into increased biomass at the cost of increase in oxidative product (here, oxygen).

1. Introduction

Life on earth, in large part, has oxygenic photosynthesis at its foundation, and much of that photosynthesis occurs in microbes. Oxygenic phototrophic microorganisms such as cyanobacteria are common in reliably lit environments, where impinging photons provide, often, a more than sufficient energy source even at low intensity, and carbon dioxide (or related chemical species) provides a reliable and abundant carbon source. When the other fundamental component of photosynthesis, water, is also available, then phototrophic based life is likely. In many cyanobacteria, nitrogen fixation can even be supported due to the abundance of photon energy. It is perhaps surprising, then, that the process of photosynthetic fixation of carbon dioxide into reduced carbon suitable for biosynthesis has, seemingly, a significant inefficiency due to the competition by oxygen for inorganic carbon binding sites, here denoted as photorespiration.

Thus, we focus on processing of inorganic carbon, i.e., carbon fixation, a central component of oxygenic phototrophy, and on its principle byproduct, molecular oxygen. Oxygenic phototrophy uses photon energy to extract electrons from water and eventually apply those electrons to fix inorganic carbon, while, in the process, oxygen is produced: electron source (effectively here, water and light) feeds electron sink (inorganic carbon) while producing oxidative byproduct (molecular oxygen). Implicit to this assembly line is the need for extracellular, macroscale transfer of inorganic carbon and of oxygen. Rates of macroscale transport (advective and/or diffusive) are largely beyond the control of individual cells and thus oxygen concentration may serve as a signal of transport limitation, triggering photorespiration. High photon flux can be an even further aggravating factor if transport of inorganic carbon into the photosynthesizing machinery cannot keep pace.

Nevertheless, the oxygenic phototrophic “business model” is generally robust and capable of being largely self sufficient. Strikingly, however, phototrophic organisms are often found in multispecies consortia together with heterotrophs. It is not immediately clear why this should be the case, as competition for resources, e.g., space or nutrients, is possible, and it seems that oxygenic phototrophs might be expected to be able to outcompete heterotrophic neighbors for those resources. Even so, multispecies communities are observed including in environments where heterotrophs might not be able to persist on their own [1]. Further, there are at least some examples of communities where resident phototrophs lack critical anabolic capabilities and must instead rely on nearby organisms to supply them [2]. Here, we explore the possible utility of interaction via organic/inorganic carbon exchange. Note there are other possible advantages in adding a heterotroph to the autotroph community. For example, heterotroph-induced oxygen usage or moderation of variation in redox potential may mitigate transport limitation.

The models presented here, both for single species (an oxygenic phototroph we call cyanobacteria) as well as for a combined two species system (cyanobacteria plus a generic heterotroph) are based in a chemostat platform. The chemostat serves as a simple and convenient way to mimic an environment where, over long times, nutrient inflow and byproduct outflow occur at rates determined by external environmental factors. From this viewpoint, a chemostat is a natural choice here due to its simplicity and also the steady oligotrophic environment it models, and thus hopefully is a reasonable bridge between abstract modeling and empirical observations.

In fact, comparison of population models with population scale observations has a well established methodology in microbial ecology. Of late, however, rapidly increasing use of molecular level technology (e.g., high throughput sequencing) has dramatically changed the nature and scale of these observations. As a result, in principle and increasingly also in practice, detailed data describing microbial capability and function is available. This information can and should potentially be used to understand how microbes exploit and alter their environment. There is a substantial gap, however, between molecular behavior at the cellular microscale and emergent community function at the population macroscale. Intermediate between the two, progress is being made in translating genomics information into models of cell dynamics [3]. Annotation of gene sequences into so-called wiring diagrams is becoming increasingly reliable and automatable. These diagrams encode cell physiology along with regulatory machinery and are accompanied by an intimidating list of unknown rate constants. However, gene encoded functions relevant to metabolic processes are naturally organized into gene pathways [4,5,6], and then, under the often reasonable assumption of steady state, balance of influx and outflux through these pathways makes choice of individual reaction rates within any particular pathway unnecessary, replaceable instead by a single flux through that entire path. Regulatory function can be characterized as a management of resource allocation between different paths and then modeled by imposing optimality criteria on that allocation [7,8]. The result is an enormous simplification: cell function is now characterized by only a limited number of rates of cellular inflow and outflow of substrates and byproducts together with an optimization principle to divide them between available metabolic pathways.

Still, there are two significant though not unrelated requirements for use of such analyses. First, despite the reduction, there remain, generally, many available and redundant metabolic flux pathways encoded by any one genome and so, as mentioned, some principle is necessary in order to decide how flux is to be distributed between those pathways. Second, also as mentioned, rates of substrate flow into the cell and byproduct flow out of the cell need to be characterized. The first of these issues couples to the second which then couples to the environment in which the cell and its community find themselves [9]. Conversely, though quantities on the large environmental scale are oftened characterized by concentrations, from this point of view it seems, rather, that fluxes are natural quantities at the cell scale. Thus, beyond the immediate aim of studying photorespiration, a further goal of this study is to suggest ways to match models at cell and population levels.

2. Materials and Methods

2.1. Model Description

We study productivity of two interacting species, one a photoautotroph with cell density ${\widehat{\mathrm{P}}}_1\left(t\right)$ and the other a heterotroph with cell density ${\widehat{\mathrm{P}}}_2\left(t\right)$, both of which are growing in a well mixed chemostat with dilution rate D [10], exposed to photon influx, see Figure 1. For simplicity, we neglect transport across the chemostat-air interface or suppose that such an interface does not exist, and let external inflow and outflow of dissolved quantities be governed by the chemostat dilution rate D. Conversely, since our aim is to study possible mutualistic or commensal effects over long times, we do not include diel light variation effects, considering them, for such purposes, to be relatively short time phenomena that can also be averaged out.

A central element of the model is the tracking of carbon flow through a microbial communtiy. As such it is convenient to measure all carbon carrying quantities in terms of carbon moles (Cmoles), e.g., to measure phototroph and heterotroph populations by the total moles of carbon they incorporate. We assume here, for convenience only, that cell sizes and densities are similar, i.e., that the total carbon moles per microbial cell, denoted c, is a constant and is the same constant for both phototrophs and heterotrophs. To convert populations from units of cells/volume to units of Cmoles/volume, we change to ${\mathrm{P}}_1\left(t\right)=c{\widehat{\mathrm{P}}}_1\left(t\right)$ and ${\mathrm{P}}_2\left(t\right)=c{\widehat{\mathrm{P}}}_2\left(t\right)$, both with units Cmoles/volume. In addition, we measure both dissolved component densities $\mathrm{IC}\left(t\right)$ (pooled inorganic carbon, Cmoles/volume), $\mathrm{OC}\left(t\right)$ (organic carbon, Cmoles/volume) in Cmoles, and ${\mathrm{O}}_2\left(t\right)$ (oxygen, molecular oxygen moles/volume) in moles of molecular oxygen. In computations, we use liters as volume units. It is assumed that oxygen concentrations always remain sufficiently low so that oxygen remains in solution and a gas phase does not occur. Note that we use the notation ${\mathrm{O}}_2$ both to denote molecular oxygen and its concentration. Inorganic carbon in solution, $\mathrm{IC}\left(t\right)$, consists of aqueous ${\mathrm{CO}}_2$ and related dissolved forms, notably aqueous bicarbonate ${\mathrm{HCO}}_3^{-}$; we do not distinguish between the forms here, though phototrophs generally do. Organic carbon here is supposed for specificity to consist of glycolic acid ${\mathrm{C}}_2{\mathrm{H}}_4{\mathrm{O}}_3$, a byproduct of photorespiration. Note that the term photorespiration has been used in the literature to designate a number of different mechanisms that, effectively, oxidize photosynthetically fixed carbon [11]. We consider here only one of those mechanisms, namely oxygenic activity of RuBisCO (ribulose bisphosphate carboxylase) secretion from the cell of partially oxidized organic carbon in the form of glycolic acid. For brevity, however, we will use the umbrella label photorespiration for this single type.

For purposes of tracking carbon, we could, as is commonly done in microbial population models, also include a microbial species ${\mathrm{Q}}_1\left(t\right)$ (Cmoles/volume) consisting of inactive photoautotroph biomaterial damaged (or killed/lysed) due to oxidative stress or in some other manner, as well as a similar heterotroph damage species ${\mathrm{Q}}_2\left(t\right)$ (Cmoles/volume). For simplicity and clarity, however, we include oxidative damage only through its direct effect on photosynthetic machinery. Note, though, that as a result, importance of oxidative damage and its amelioration are, if anything, likely underestimated in the later results.

The general form of the equations used here for a chemostat with photon flux ν are

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{P}}_1& =& (\eta g_1(\mathrm{IC},{\mathrm{O}}_2;\nu )-D){\mathrm{P}}_1,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{P}}_2& =& (g_2(\mathrm{OC},{\mathrm{O}}_2)-D){\mathrm{P}}_2,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{IC}& =& -Y_{{\mathrm{P}}_1,\mathrm{IC}}^{-1}{g}_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1+y_{{\mathrm{P}}_2,\mathrm{IC}}^{-1}{g}_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2+D({\mathrm{IC}}_0-\mathrm{IC}),\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{OC}& =& Y_{\mathrm{IC},\mathrm{OC}}^{-1}(1-\eta )g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1-y_{{\mathrm{P}}_2,\mathrm{OC}}^{-1}{g}_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2-D\mathrm{OC},\hfill \\ \hfill \frac{d}{dt}{\mathrm{O}}_2& =& \left((Y_{\mathrm{IC},{\mathrm{O}}_2}^{-1}-Y_{\mathrm{OC},{\mathrm{O}}_2}^{-1}(1-\eta ))g_1(\mathrm{IC},{\mathrm{O}}_2;\nu )\right){\mathrm{P}}_1\hfill \end{array}$$

(4)

$$\begin{array}{ccc}& & -y_{{\mathrm{P}}_2,{\mathrm{O}}_2}^{-1}{g}_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2+D(O_{2,0}-{\mathrm{O}}_2),\hfill \end{array}$$

(5)

where the various subscripted $Y_{\alpha ,\beta }$’s (associated with P${}_1$) and $y_{\alpha ,\beta }$’s (associated with P${}_2$) are yield coefficients, all of which are fixed by stoichiometry, with units of Cmoles of α per Cmoles of β or moles of ${\mathrm{O}}_2$. The parameters $k_1$ and $k_2$ indicate specific rates of deactivation of active biomaterial and could be functions of O${}_2$. The function $\eta =\eta (\mathrm{IC},{\mathrm{O}}_2)$ is related to photorespiration, and will be defined later. Terms containing rate $g_1$ are involved in the photobiosynthesis and/or photorespiration pathways and terms containing rate $g_2$ are involved in the heterotrophic biosynthesis pathway. All internal metabolic rates are fixed by the three pathway (phototroph biosynthesis, photorespiration, and heterotroph biosynthesis) rates so that they need not be parameterized in detail except through the single rate functions $g_1$ and $g_2$ together with branching parameter η: this is a consequence of the powerful assumption of short timescale equilibration of metabolic pathways [6]. For easy reference, see

Table 1. State variables (left) and key environmental parameters (right) for system (1)–(5).

State Quantities			Key Environmental Parameters
Symbol	Description	Units	Symbol	Description	Units
P${}_1$	Phototroph Concentration	Cmol·L${}^{-1}$	D	Dilution Rate	s${}^{-1}$
P${}_2$	Heterotroph Concentration	Cmol·L${}^{-1}$	ν	Photon Flux	$\mathsf{\mu }$E m${}^{-2}$·s${}^{-1}$
IC	Inorganic Carbon Concentration	Cmol·L${}^{-1}$	IC${}_0$	Inflow IC Conc.	Cmol·L${}^{-1}$
OC	Organic Carbon Concentration	Cmol·L${}^{-1}$
O${}_2$	Oxygen Concentration	Omol·L${}^{-1}$	O${}_{2,0}$	Inflow O${}_2$ Conc.	Omol·L${}^{-1}$

and

Table 2. State variables (left) and key environmental parameters (right) for system (1)–(5).

Rate Functions			Yield Parameters
Symbol	Description	Units	Symbol	Description	Reference
η	Photorespiration Branching Function	–
$g_1$	Photosynthesis Rate	s${}^{-1}$	$Y_{\alpha ,\beta }$	Phototroph Yields α per β	see Section 2.3.1
	Photobiosynthesis Rate $=\eta g_1$
	Photorespiration rate $=(1-\eta )g_1$
$g_2$	Heterotroph Biosynthesis Rate	$s^{-1}$	$y_{\alpha ,\beta }$	Heterotroph Yields α per β	see Section 2.4

. Details for individual terms in (1)–(5) will be provided below.

The system energy is supplied through the photon flux ν which serves to drive carbon reduction through photosynthesis. We assume that microbe populations are sufficiently sparse in the chemostat so that no significant shading occurs, though one could introduce a shaded photon flux ${\nu }_{\mathrm{shade}}={\nu }_{\mathrm{shade}}(\nu ,{\mathrm{P}}_1,{\mathrm{P}}_2,{\mathrm{Q}}_1,{\mathrm{Q}}_2)$ (note that even in the case of shading, all microbes effectively receive the same photon flux over time due to the well mixed assumption). The other environmental conditions included are the dissolved concentrations IC, OC, and ${\mathrm{O}}_2$. Inorganic carbon and oxygen flow into the chemostat at concentrations ${\mathrm{IC}}_0$ and $O_{2,0}$, respectively. Exchange of O${}_2$ and CO${}_2$ with the atmosphere is neglected but inclusion would not be expected to change qualitative conclusions. The inflow is assumed to be free of organic carbon. Note that non-negative initial conditions are required for all quantities but have only transient influence, on a $D^{-1}$ time scale, except/unless ${\mathrm{P}}_1\left(0\right)=0$ (in which case ${\mathrm{P}}_1\left(t\right)=0$ for all t) or ${\mathrm{P}}_2\left(0\right)=0$ (in which case ${\mathrm{P}}_2\left(t\right)=0$ for all t). Hence, later, we will ignore transients and study steady states.

Photosynthesis drives ecology through conversion of photons to chemical energy (photons power $\mathrm{ADP}\to \mathrm{ATP}$, say) but also, and possibly more importantly, through production of reducing power, referred to here as electrons. In fact, we will not consider energy production and, rather, implicitly track electrons through degree of reduction (see Appendix A) as the more important quantity. A key step in oxygenic photosynthesis is the splitting of ${\mathrm{H}}_2\mathrm{O}$ into, for our purposes, a combination of ${\mathrm{O}}_2$ and reducing power. Oxygen’s importance goes beyond its role as reactant; it also is an important contributor to degree of reduction balance of the entire oxygenic photosynthesizing system. In fact, in the model presented here, oxygen is the only explicit quantity with negative degree of reduction and hence, by proxy, its concentration is central to community redox state and hence to community function.

2.2. Metabolic Pathways

From an engineering point of view, organism metabolics operate somewhat like chemical processing networks so that they and implicitly resulting ecological interactions, are conveniently represented in terms of what are called metabolic pathways, chains of reactions that convert external substrates into external byproducts (though cycles of internal reactions might also be considered as pathways). Organisms themselves might be viewed as collections of such reaction chains, interacting with each other while producing fluxes at rates which must be consistent with external flux constraints. For example, in the case of a simple chemostat, external inflow and outflow fluxes are set by dilution rate D. While we look here to adopt the point of view of organism metabolisms as collections of pathways, at the same time we want a simple system able to illustrate basic principles of a phototroph-heterotroph interactions. Thus, while detailed metabolic models exist including for cyanobacteria [12,13,14,15] as well as for communities [16,17], we reduce system metabolics to the interaction of three particular pathways: photosynthesis-driven biosynthesis, photorespiration (in a restricted sense as previously noted) in the phototrophs, and aerobic respiration-driven heterotrophic biosynthesis. Community function is determined by the rates at which these pathways operate; the environment, through chemostat inflow and outflow, constrains community function by constraining these rates, though the community, specifically here the photoautotrophic cyanobacteria, have some freedom to choose them.

There are two specific rate functions in system (1)–(5): rate of carbon fixation $g_1(\mathrm{IC},{\mathrm{O}}_2;\nu )$ of the photoautotrophs via the photosynthesis/photorespiration pathways and rate of growth $g_2(\mathrm{IC},{\mathrm{O}}_2)$ of heterotrophs via a pathway for catabolysis of available organic carbon. In addition, there is a branching parameter η which determines percent of fixed carbon going to photoautotroph growth versus photorespiration. Between the three $g_1$, $g_2$, and η, the rates of the three pathways are determined. In balance, there are essentially two types of constraints: (1) cellular inflow rates of photons and inorganic carbon, with photosynthesis determined by the minimum rate of the two, as well as (2) system inflows and outflows determined by the chemostat dilution rate D. Cells must be able to synchronize system rates, dilution rate here, to direct pathway inputs, photon and inorganic carbon in the case of photosynthesis, and oxygen and organic carbon in the case of heterotrophic anabolysis. In principle, cellular outflow rates for pathway products may also constrain, but we suppose here, for the particular pathways studied, that these rates are essentially free.

Relation to Pathway Analysis

Metabolic network analysis of a system of m metabolites with internal concentrations $c_i$, $1\le i\le m$, and n reactions with rates $v_j(c_1,c_2,...,c_m)$, $1\le j\le n$, starts from a metabolic map that can be represented by a set of equations of the form

$$\frac{d{c}_i}{dt}=\sum _j{N}_{ij}{v}_j$$

where $N_{ij}$ is a stoichiometric coefficient, possibly negative, for production of metabolite i via reaction j. A rate $v_j$ can be determined as a function of the concentrations $c_i$ and is parameterized by rate constants. These rate constants are often unknown, but if steady state is assumed then the problem reduces to characterization of the null space

$$N\mathbf{v}=\mathbf{0}$$

(6)

of the $m\times n$ stoichiometric matrix N in a useful way by somehow identifying important pathway vectors $\mathbf{v}$ from this null space. (Precisely, a pathway consists of the reactions corresponding to non-zero entries in a pathway vector $\mathbf{v}$; a pathway vector encodes the flux through each of those reactions). Note that knowledge of rate constants is unnecessary to solve the steady state Equation (6) and thus also unnecessary to determine pathway vectors, though steady state internal concentrations $c_i$ cannot be computed without these rate constants.

One objective here is to proceed a further step by connecting internal metabolic activity, as encoded by those distinguished pathways, to community dynamics, e.g., connecting information extracted from (6) to the community model, as stated in (1)–(5). To do so, we use the (significantly) reduced metabolic maps as shown in Figure 2, explained in detail later. The interiors of the dashed domains in Figure 2 correspond to the interiors of the circled objects P${}_1$ and P${}_2$ of Figure 1. Circled objects in Figure 2 are generalized metabolites and arrowed curves are generalized reactions. Metabolites that are associated with reactions exiting a dashed domain are “seen” by the environment and hence explicitly tracked in the model (1)–(5); other metabolites are internal (in this case, only electrons e) so not directly observed in the environment and thus not explicit in the model, i.e., do not have tracking equations. Interior dynamics are assumed to be at quasi-steady state, that is, are able to quickly equilibrate to time on the community interaction time scale. Later, we will suppose a third, longer time scale on which the community also reaches steady state.

We do not construct here the stoichiometric matrix corresponding to the metabolic network in Figure 2, but rather proceed directly to its elementary flux modes [6] which mathematically are, where reversible reactions are not present as is the case here, non-negative solutions of (6) for which no other non-negative solution containing a proper subset of non-zero entries exists. That is, an elementary flux mode is, roughly, a realizable pathway through the metabolic network that does not contain within itself any smaller realizable pathways. Non-negative linear combinations of the complete list of all elementary flux modes of a given network generate all allowable solutions of (6). For the system in Figure 2, there are three elementary modes, see Figure 3, corresponding to (1) biomass production and (2) photorespiration in the phototrophs and to (3) biomass production in the heterotrophs. All realizable steady states of the system can be uniquely written as positive combinations of these three elementary modes.

The metabolic state of the two population of microbes is thus described by the rates of the three elementary modes. These rates are determined by the flux rates into cells of mode inputs and out of cells of mode outputs, all of which generally depend on external concentrations of those inputs and outputs. Note that, for a given mode, setting any one of the input or output rates determines all reaction rates through the entire mode due to stoichiometric constraints and the steady state assumption. Hence, though many individual reactions are involved, each with its own rate constant, the three modes can be described by only three rates in total (one for each).

It should be noted that, for the sake of simplicity, we in fact compromise flux mode balance in one respect: flux balance through the phototroph electron compartment (circled e in Figure 2. left) is not explicitly enforced in the carbon limited case when the supply of inorganic carbon is insufficient to match the electron supply. Instead, in that particular case, some excess electrons are removed from the system. However, we keep track of their flux implicitly through photoinhibition – those excess electrons effectively combine with biomass and oxygen to remove some biomass. To maintain explicit flux balance, we would need to track them say to reaction with reactive oxygen species or wherever else they may go. As we implicitly suppose that such products leave the system, explicit tracking would complicate the model without advantage, as particular mechanisms of excess electron removal and damage are not the principle focus here.

The two phototroph modes operate in parallel and hence compete directly, in a sense, for inputs (inorganic carbon and photons). They operate in series with the single heterotroph mode and interact with it indirectly through external concentrations of dissolved quantities. In the case of a chemostat with dilution rate D, transport in and out of the chemostat of all quantities also proceeds at rate D. It must thus necessarily be the case that biosynthesis modes also operate at rate D (or else at rate 0) placing two constraints on mode rates. Hence we have one remaining condition determined, that of photorespiration. From there, concentrations of external quantities are determined by consistency with mode rates. These external concentrations effectively determine steady state biomass concentrations.

2.3. Photosynthesis

Microbial oxygenic photosynthesis can be divided into two steps, the light reaction followed by the dark reaction (also called the Calvin cycle), so named because photons are involved only in the first step [11]. The entire process uses energy from incoming photons to split ${\mathrm{H}}_2\mathrm{O}$ producing ${\mathrm{O}}_2$ and electrons, which, in the form $\mathrm{NADPH}$, are used to fix ${\mathrm{CO}}_2$. The light reaction, which is the oxygenic step, can be summarized for our purposes by

$$2{\mathrm{H}}_2\mathrm{O}+8\mathrm{photons}⟶4{\mathrm{H}}^{+}+4\mathrm{electrons}+{\mathrm{O}}_2$$

(7)

and the dark reaction, the carbon fixation step, can be summarized, again for our purposes, by

$${\mathrm{CO}}_2+\frac{1}{2}(1-\eta ){\mathrm{O}}_2+\omega \mathrm{electrons}⟶\eta {\mathrm{CH}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}+(1-\eta ){\mathrm{CH}}_2{\mathrm{O}}_{1.5},$$

(8)

with both formulas balanced for carbon and degree of reduction (and only carbon and degree of reduction, for simplicity). Note that the dark reaction also consumes energy in the form of ATP, which may also be of importance for cellular energy balances, but ATP cycling is not considered here. Model parameters associated with (7)–(8) are described in

Table 3. Key photosynthesis-related functions and parameters.

Symbol	Description	Units	Definition
A	Average cell cross-sectional area	$\mathsf{\mu }$m${}^2$	Appendix B
c	Carbon moles per cell	Cmole/cell	Appendix B
e	Electron production rate by the light reaction	emole/Cmole·s	Equation (12)
f	Net oxygen per photosynthetically fixed carbon	Omole/Cmole	Equation (10)
I	Photoinhibition function	–	Equation (15)
α	photosynthesis efficiency factor	–	Equation (12)
${\gamma }_1$	inverse specificity factor	Cmole/Omole	Equation (16)
${\gamma }_2$	excess electron capacity	s/Omole	Equation (15)
ϵ	Maximum electron consumption rate	1/s	Equation (13)
η	RuBisCO inorganic carbon binding probability	–	Equation (16)
ν	environmental photon flux	$\mathsf{\mu }$E/m${}^2·$s	–
ω	Electron demand: emoles needed to fix a cmole	emole/Cmole	Equation (9)

.

The division into two steps as formulated by (7) and (8) has important consequences. In particular, at least within the formulation of the model, the light reaction is governed by photon supply (H${}_2$O being assumed to be abundant) whereas the dark reaction rate is determined by both output rate of the light reaction (in electrons) as well in inflow rate of inorganic carbon. Hence, effectively, the carbon fixation rate function $g_1$ is determined by the rate at which the dark reaction proceeds, which may be slower than that of the light reaction in the case of limiting CO${}_2$ (or, for us, IC). If so, the excess electrons effectively recombine with oxygen or some other oxidant, but may cause inhibition of photosynthesis in the process.

A second point of note is that photosynthesis as formulated by (7) and (8) is divided into two separate pathways weighted by the branching function η, with $0\le \eta \le 1$. The valuation $\eta =1$ corresponds to all fixed carbon being used for a growth pathway, while $\eta <1$ indicates that some of the fixed carbon is instead allocated to a pathway that results in secretion of organic carbon from the cell (with $\eta =0$ corresponding to all fixed carbon going to the secretion pathway, though this outcome would not allow a viable population). Within the model presented here, there are seemingly two apparent advantages to $\eta <1$ over $\eta =1$: first, $\eta <1$ results in consumption of oxygen, see (8), which may alleviate effects of oxygen oversupply, and second, $\eta <1$ results in secreted organic carbon which can be used to supply a population of heterotrophs which in turn produce inorganic carbon as a byproduct. These points are discussed in more detail below, but note that we suppose here that the value of η is subject to control by cyanobacterial cells themselves. Hence it may be that η is in some way regulated through an optimization process.

2.3.1. Carbon Fixation

The righthand side of (8) consists of two types of fixed organic carbon, each produced as a consequence of Calvin cycle reactions, a carboxylase reaction and an oxygenase reaction. That is, formula (8) actually combines contributions from two pathways: production of ${\mathrm{CH}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}$ via biosynthesis, weighted by η, and production of ${\mathrm{CH}}_2{\mathrm{O}}_{1.5}$ via photorespiration, weighted by $1-\eta$.

New biomass is approximated as ${\mathrm{CH}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}$ [18]. In actuality, the dark reaction only produces a precursor (glyceraldehyde 3-phosphate) and biosynthesis is completed elsewhere, but for purposes of electron balance, it is convenient to use the biomass proxy formula ${\mathrm{CH}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}$ as the ultimate biomass output, see the source term in (1). ${\mathrm{CH}}_2{\mathrm{O}}_{1.5}$ (carbon-normalized glycolate) is a soluble byproduct of photorespiration and is assumed to be excreted from the cell; we track its concentration in the chemostat as the quantity $\mathrm{OC}\left(t\right)$, see (4). Recall that, for purposes of representation in terms of carbon moles, all carbon compounds are normalized so that the number of carbon atoms is one.

Given the branching function η, described and parameterized below, then

$$\omega =\omega \left(\eta \right)=4.7\eta +5(1-\eta )=5-0.3\eta$$

(9)

is the electron demand (emole/Cmole), the number of moles of electrons needed to fix a mole of inorganic carbon (ω is related to carbon-oxygen-demand, a quantity sometimes used in engineering applications). We use degree of reduction 4.7 for ${\mathrm{CH}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}$ (autotroph) and degree of reduction 3.0 for ${\mathrm{CH}}_2{\mathrm{O}}_{1.5}$, see Appendix A. The coefficient in (9) of $(1-\eta )$ is 5 rather than 3, though 3 is the degree of reduction of ${\mathrm{CH}}_2{\mathrm{O}}_{1.5}$, because the left-hand side term $(1/2)(1-\eta ){\mathrm{O}}_2$ in (8), with degree of reduction $-2(1-\eta )$, effectively transfers to the righthand side for purposes of computing electron demand. For reference, note that the yield coefficients of moles of molecular oxygen per fixed Cmole are computed from (7) and (8) to be

$$\begin{array}{ccc}\hfill Y_{\mathrm{IC},{\mathrm{O}}_2}^{-1}& =& \frac{\omega \left(\eta \right)}{4},\hfill \\ \hfill Y_{\mathrm{OC},{\mathrm{O}}_2}^{-1}& =& \frac{1}{2},\hfill \end{array}$$

with rates proportional to $g_1{\mathrm{P}}_1$, see the righthand side of (5). The coefficient of the first term in the righthand side of (5),

$$f\left(\eta \right)=Y_{\mathrm{IC},{\mathrm{O}}_2}^{-1}-Y_{\mathrm{OC},{\mathrm{O}}_2}^{-1}(1-\eta )=\frac{3}{4}+\frac{17}{40}\eta ,$$

(10)

indicates the net yield of oxygen moles per mole of photosynthetically fixed carbon and is important in the results presented here. Note that by varying η between 0 (all fixed carbon goes to excreted, soluble carbon) and 1 (all fixed carbon goes to new biomass), organisms vary oxygen production by about 60%. Of course, this should be understood as at best a rough estimate, since the presented model greatly simplifies the true biochemistry, but nevertheless the variation is potentially significant. At the same time, electron demand is much less sensitive, varying only by about 6% as η varies between 0 and 1, so it might seem that photorespiration provides a means to reduce oxygen production without significantly reducing capacity to process the photosynthetically driven electron stream.

Electrons are not normally free in solution but, rather, are transported by carrier compounds, particularly ${\mathrm{NADP}}^{+}/\mathrm{NADPH}$, and passed along through redox reactions [19]. As a convenience, though, we account for electrons directly rather than track ${\mathrm{NADP}}^{+}$ and $\mathrm{NADPH}$. Note that in each reaction Cmoles and degree of reduction (see Appendix A) are balanced as the key governing quantities. Other reaction components, e.g., N, are considered of secondary importance and are not balanced. Doing so would require introduction of more reactions, obscuring the main points. For example, protons are in excess and can be assumed to be buffered by the aqueous environment through mechanisms not of direct importance to the modeling aims here. Note in passing, though, that in some instances growth is limited by availability of quantities other than those tracked here, e.g., by limitation in fixed nitrogen. We suppose that this is not the case here.

Reactions (7) and (8) together comprise the photosynthesis and photorespiration pathways (which branch from each other in the dark reaction step), with ${\mathrm{CH}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}$ being the output of the photobiosynthesis pathway and ${\mathrm{CH}}_2{\mathrm{O}}_{1.5}$ the output of the photorespiration pathway. Excess molecular oxygen is also an output of both. An emphasis on rate rather than concentration is key and all internal reaction rates are effectively slaved to rates of inflow and outflow to/from the cell. The only other needed parameters are the stoichiometric ones, which are known from the pathway descriptions, in this case (7) and (8). Hence it is important to characterize governing rates, particularly those that have limiting or other important roles. The principle inputs of interest to photosynthesis are photons and ${\mathrm{CO}}_2$ (water is plentiful at least in a chemostat) and growth rate is limited by the lesser availability of the two. We assume that the principle bottleneck for ${\mathrm{CO}}_2$ inflow is transport (more specifically, here, transport of inorganic carbon – recall that we do not distinguish between inorganic carbon species) from outside the cell to the photosynthetic machinery inside. As is commonly practice, e.g., [20], we approximate this transport rate by a Michaelis-Menten function of the form

$$\delta \left(\mathrm{IC}\right)=r_{\mathrm{IC}}\frac{\mathrm{IC}}{K_{\mathrm{IC}}+\mathrm{IC}}\left({\mathrm{time}}^{-1}\right)$$

(11)

where $r_{\mathrm{IC}}$ and $K_{\mathrm{IC}}$ are, respectively, maximal transport rate and half-saturation of cross-membrane transport, see Appendix B. There are a number of mechanisms cells can use to influence transport, notably carbon capturing and active transporters [11,21]. From our point of view, carbon capturing effects can, roughly, be replaced by decrease in the inverse specificity factor ${\gamma }_1$ defined below, and active transporters can influence parameters $r_{\mathrm{IC}}$ and $K_{\mathrm{IC}}$. Ultimately transport rate of inorganic carbon into cells over time is limited by concentrations outside of the cell and, more particularly, transport rates of inorganic carbon into and out of the local environment.

Rate δ then needs to be compared to the rate at which the dark reaction (8) can use the electrons to match with the inflowing inorganic carbon. The light reaction (7) provides that electron supply. Photons flow through the chemostat with constant flux ν (photons/area·time) set externally as a parameter, but in actuality enter photosynthesis machinery at an effective rate ${\nu }_{\mathrm{eff}}=A\alpha \nu$ (photons/cell·time) where A is cell cross-section (area/cell) and $0\le \alpha \le 1$ is an efficiency factor (unitless), see Appendix B. The parameter alpha accounts for photons that impact the cell but do not result in oxygen and electron production, either because they do not enter the photosynthetic process at all or because their end impact is shunted to non-photochemical quenching precesses such as Mehler reactions. These latter mechanisms may have other outcomes such as ATP generation which are supposed here to be non-limiting (though can have negative impacts at high enough levels) so are not considered. Note that alpha can be scaled into nu, so changing efficiency is equivalent, in the model, to changing light intensity. Note that α in fact measures of the efficiency of the process of electrons impacting the cell all the way to production of electrons and reductant, and in principle may be a function of conditions such as photon flux, oxygen concentration, etc., though we do not try to model these effects here. The light reaction component of photosynthesis then produces electrons at rate ${\mathrm{Y}}_{\nu }{\nu }_{\mathrm{eff}}$ (electrons/cell·time), with yield ${\mathrm{Y}}_{\nu }=1/2$ electron/photon. Electron production rate by the light reaction, per Cmole of biomass, is thus

$$e=\frac{Y_{\nu }{\nu }_{\mathrm{eff}}}{c}=\frac{A\alpha \nu }{2c}(\mathrm{emole}/\mathrm{Cmole}·\mathrm{time})$$

(12)

(recall that c is the number of Cmoles per cell). The dark reaction then processes these electrons together with inorganic carbon. Note however, from (8) and (9), that this processing depends on the division of output between biomass and soluble inorganic production. Hence the rate at which electrons are actually consumed by the dark reaction (at least if inorganic carbon supply is not limited) is

$$ϵ=\frac{e}{\omega }\left({\mathrm{time}}^{-1}\right),$$

(13)

which can be understood as the maximum rate at which photon flux can drive carbon fixation through combined photobiosynthesis and photorespiration.

In fact, reduction of inorganic carbon proceeds by matching, in a sense, incoming inorganic carbon with incoming electrons, with rates set by $\delta \left(\mathrm{IC}\right)$ and ϵ, respectively. Since $\delta \left(\mathrm{IC}\right)\ne ϵ$ in general, then in fact reduction can proceed at best at rate $min(\delta ,ϵ)$. In the spirit of rate-based modeling, we suppose this minimum to largely govern the actual reduction rate, so that photosynthesis rate, more particularly, the dark reaction rate, is

$$g_1(\mathrm{IC},{\mathrm{O}}_2;\nu )=\left\{\begin{array}{ccc}ϵ\hfill & ϵ-\delta \le 0\hfill & \left(\mathrm{light}\mathrm{limited}\right)\hfill \\ \delta I\hfill & ϵ-\delta >0\hfill & \left(\mathrm{carbon}\mathrm{limited}\right)\hfill \end{array}\right.$$

(14)

where $I=I(ϵ-\delta ,{\mathrm{O}}_2)$ is a photoinhibition function, defined below, of excess electrons should there be any. Note that O${}_2$ dependence in $g_1$ arises from O${}_2$ dependence in ω and I.

2.3.2. Photoinhibition and Oxidative Stress

In the case that $ϵ>\delta$, i.e., the rate of the normalized electron production is greater than the rate of inorganic carbon inflow, excess electron production can lead to inhibition of photosynthesis machinery and other apparatus via saturation of electron transport structure and consequent formation of harmful radical oxygen species as well as other undesirable effects [11,22]. These effects have been modeled with an inhibition function [23,24] which allows for removal of excess electrons without detriment, to a point, after which reduction in growth rate occurs [22,25]. These inhibition models, however, are generally functions only of photon flux rate and not, for example, dependent on IC and O${}_2$ concentrations or transport rates, though such dependence is likely important, at very least through any mismatch of electron production rate with inflow rate of electron donors, and is certainly central to the model presented here. Thus, we define a photoinhibition function to take account of electron-IC mismatch of the form

$$I(ϵ-\delta ,{\mathrm{O}}_2)=\left\{\begin{array}{cc}1\hfill & ϵ-\delta \le 0\hfill \\ {(1+{\gamma }_2^{-2}{(ϵ-\delta )}^2{\mathrm{O}}_2^2)}^{-1}\hfill & ϵ-\delta >0\hfill \end{array}\right.$$

(15)

where ${\gamma }_2$ is an excess capacity coefficient. (Note, from (14) that only the $ϵ-\delta >0$ definition is relevant.) The quadratic dependencies on O${}_2$ and $ϵ-\delta$ are ad hoc forms meant to model the capacity for cells to avoid or repair damage of small mismatches; small values of $ϵ-\delta$ and O${}_2$ do not result in significant net damage whereas large values might. The parameter ${\gamma }_2$ is chosen to provide a reasonable high light/oxygen cutoff on growth. As a consequence of the degree of arbitrariness in form and parameterization of I, we avoid conclusions which would appear to rely on its particulars beyond a general tendency to inhibit growth in carbon limited conditions.

We note here that at moderately high concentrations O${}_2\gtrsim 5·{10}^{-4}$ Omoles/L, oxygen may come out of solution, providing effectively a method for limiting effects of oxygen stress. Such critical oxygen concentrations are not reached in computations shown here, but can occur at environmentally reasonable light intensities in some situations.

2.3.3. Photorespiration

A key step in dark reaction carbon fixation is binding of ${\mathrm{CO}}_2$ to the enzyme RuBisCO. However, as it happens, ${\mathrm{O}}_2$ competes for the same binding site as ${\mathrm{CO}}_2$, and when a molecule of ${\mathrm{O}}_2$ does in fact bind then glycolate ($\mathrm{C}{\mathrm{H}}_2{\mathrm{O}}_{1.5}$, degree of reduction $+3$) is produced in the stead of further reduced biomaterial ($\mathrm{C}{\mathrm{H}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}$, degree of reduction $+4.7$ for phototrophs). We refer to this process as photorespiration (though as noted earlier, photorespiration can be used as an umbrella term for a number of re-oxidixing processes). We denote the probability of ${\mathrm{CO}}_2$ binding to RuBisCO by η, with

$$\eta =\frac{a_c\mathrm{IC}}{a_c\mathrm{IC}+a_o{\mathrm{O}}_2}=\frac{1}{1+{\gamma }_1\frac{{\mathrm{O}}_2}{\mathrm{IC}}}$$

(16)

where $a_c$, $a_o$ are binding affinities and ${\gamma }_1$, the ratio of those affinities, is the inverse specificity factor (with respect to IC versus O${}_2$). Recall that we confuse inorganic carbon concentration IC here with CO${}_2$ concentration, supposing that inorganic carbon in forms other than CO${}_2$ can readily be converted into CO${}_2$ via carbonic anhydrase enzymes.

Effectively, η is a branching function of O${}_2$ and IC that determines how much photosynthetic product goes to synthesis of new biomaterial and how much to synthesis of soluble, excretable, organic carbon. Phototrophs may have a degree of control over the value of η either directly through the structure of RuBisCO itself [26,27] or through indirect machinery such as carbon capture mechanisms, so we treat ${\gamma }_1$ as a tunable parameter and study effects of its variation.

The purpose of photorespiration (oxygenase activity of RuBisCO, to be precise), if there is one, is uncertain. It is sometimes argued to be wasteful, e.g., [20], and possibly a relic of early earth history when levels of ${\mathrm{CO}}_2$ were much higher than today, and levels of O${}_2$ lower, so that the the ratio O${}_2$/CO${}_2$ was presumably small. However, observations suggest it is not superfluous [28] and the orders of magnitude variability of γ across different species [11,26,27] suggests that there may be selective pressure at work. Photorespiration diverts carbon fixing power away from new biomass, but also note in fact the following: though glycolate has a lower degree of reduction ($+3$) than biomaterial ($+4.7$), its production requires 1/2 ${\mathrm{O}}_2$ mole per Cmole of glycolate and hence, balancing electrons, also removes an additional two electrons per glycolate. Thus, effectively, each Cmole of glycolate produced removes 5 electron moles from the system, more than the 4.7 electron moles removed per Cmole of biomaterial produced. Thus photorespiration serves to reduce electron pressure, particularly when oxygen pressure is high. At the same time, oxygen pressure is reduced. Also, photorespiration produces a supply of dissolved, reduced organic carbon, allowing the possibility of supplying a heterotroph population. Hence, accidental or not, photorespiration may have significant effects on population dynamics.

2.3.4. Fixation Stability

As a technical point that will be repeatably useful below but also seems reasonable biologically, we impose the condition

$${\eta }_{\mathrm{IC}}{g}_{1,{\mathrm{O}}_2}-{\eta }_{{\mathrm{O}}_2}{g}_{1,\mathrm{IC}}\ge 0$$

(17)

(with subscripts IC and O${}_2$ denoting partial derivatives with respect to those quantities)) which, with η as in (16), reduces to $\mathrm{IC}{g}_{1,\mathrm{IC}}\ge -{\mathrm{O}}_2{g}_{1,{\mathrm{O}}_2}$. (In fact we will only really require (17) to hold at steady state). This condition can be appreciated through linearization of the carbon fixation process applied to inorganic carbon and oxygen, i.e., linearization of the subsystem

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{IC}& =& -Y_{{\mathrm{P}}_1,\mathrm{IC}}^{-1}{g}_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1\hfill \\ \hfill \frac{d}{dt}{\mathrm{O}}_2& =& f\left(\eta \right)g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1\hfill \end{array}$$

around a state $\left({\mathrm{P}}_1,\mathrm{IC},{\mathrm{O}}_2\right)$, with associated Jacobian matrix

$$J=\left(\begin{array}{cc}\hfill -Y_{{\mathrm{P}}_1,\mathrm{IC}}^{-1}{g}_{1,\mathrm{IC}}{\mathrm{P}}_1& \hfill -Y_{{\mathrm{P}}_1,\mathrm{IC}}^{-1}{g}_{1,{\mathrm{O}}_2}{\mathrm{P}}_1\\ \hfill {\left(f{g}_1\right)}_{\mathrm{IC}}{\mathrm{P}}_1& \hfill {\left(f{g}_1\right)}_{{\mathrm{O}}_2}{\mathrm{P}}_1\end{array}\right).$$

The eigenvalues of J have non-positive real part as long as derivatives with respect to IC are non-negative, derivatives with respect to O${}_2$ are non-positive, and condition (17) holds. In the case that (17) is false, then J has an unstable eigendirection that corresponds to an instability in the fixation process: a simultaneous increase in IC and O${}_2$ levels can lead to simultaneous decrease in net fixation rate and in photorespiration, thus further amplifying IC and O${}_2$ levels, etc. Such dynamics are unsustainable. Equivalently, it can be seen that, if (17) is false, then an increment in available inorganic carbon actually reduces photosynthesis rate, see Appendix C.

Condition (17) is satisfied for reasonable choices of η and $g_1$, with one caveat, see below. We divide into two cases based on (14). In the light limited regime,

$${\eta }_{\mathrm{IC}}{g}_{1,{\mathrm{O}}_2}-{\eta }_{{\mathrm{O}}_2}{g}_{1,\mathrm{IC}}={\eta }_{\mathrm{IC}}{ϵ}_{{\mathrm{O}}_2}-{\eta }_{{\mathrm{O}}_2}{ϵ}_{\mathrm{IC}}=0,$$

satisfying (17). In the carbon limited regime

$${\eta }_{\mathrm{IC}}{g}_{1,{\mathrm{O}}_2}-{\eta }_{{\mathrm{O}}_2}{g}_{1,\mathrm{IC}}={\eta }_{\mathrm{IC}}{\left(\delta I\right)}_{{\mathrm{O}}_2}-{\eta }_{{\mathrm{O}}_2}{\left(\delta I\right)}_{\mathrm{IC}}={\eta }_{\mathrm{IC}}\delta I_{{\mathrm{O}}_2}-{\eta }_{{\mathrm{O}}_2}{\left(\delta I\right)}_{\mathrm{IC}}$$

(18)

The first term on the far right hand side is generically non-positive, while the second is generically non-negative. Note the key controlling function, $I_{{\mathrm{O}}_2}$, indicates the rate at which increasing oxygen levels increases oxidative stress; only if this rate is too large can (18) be negative. Otherwise, fixation stability condition (17) also holds in the carbon limited regime.

2.4. Heterotrophic Biosynthesis

The third pathway in the model system is a simplified heterotrophic anabolysis described by

$$2{\mathrm{CH}}_2{\mathrm{O}}_{1.5}+0.475{\mathrm{O}}_2⟶{\mathrm{CH}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}+{\mathrm{CO}}_2,$$

(19)

with stoichiometry constrained to balance carbon and degree of reduction (using degree of reduction of $\left(\mathrm{C}{\mathrm{H}}_{1.7}{\mathrm{O}}_{0.5}{\mathrm{N}}_{0.2}\right)=+4.1$ for heterotrophs, see Appendix A). As with the photosynthesis pathway, oxygen and nitrogen are not balanced; to do so would require introduction to the model of new details of secondary interest. Note that the stoichiometry determines yield coefficients $y_{{\mathrm{P}}_2,\mathrm{OC}}=1/2$, $y_{{\mathrm{P}}_2,{\mathrm{O}}_2}=4/1.9$, and $y_{{\mathrm{P}}_2,\mathrm{IC}}=1$ in Equations (3)–(5).

Reaction (19) indicates that organic carbon in the form of glycolate is further reduced to biomaterial. The increased degree of reduction is accomplished by sacrificing some of the glycolate for its electrons, a portion of which go to biomaterial and a portion of which are shunted off to carbon dioxide to maintain carbon balance.

The rate at which (19) proceeds is given by $g_2\left({\mathrm{OC},\mathrm{O}}_2\right)$ as

$$g_2\left(\mathrm{OC}\right)=r_{h}min\left(\frac{\mathrm{OC}}{K_{\mathrm{OC}}+\mathrm{OC}},\frac{{\mathrm{O}}_2}{K_{{\mathrm{O}}_2}+{\mathrm{O}}_2}\right),$$

(20)

based on the assumption that the rate of biosynthesis is controlled by the minimum rate at which biosynthesis components can be transported to biosynthesis machinery. Note as well that as a result of reaction (19), a source term for P${}_2$ appears in (2) and sink terms for OC and O${}_2$, appear in (4) and (5).

2.5. Equations

Pathway stoichiometry can now be incorporated into Equations (1)–(5). Yields $Y_{\alpha ,\beta }$ parameterize autotroph pathways (and damage) and yields $y_{\alpha ,\beta }$ parameterize the heterotroph pathway (and damage). In units of Cmoles and oxygen moles, $Y_{{\mathrm{P}}_1,{\mathrm{Q}}_1}=y_{{\mathrm{P}}_2,{\mathrm{Q}}_2}=1$ (see Section 2.3.2), $Y_{{\mathrm{P}}_1,\mathrm{IC}}=1$ (see Section 2.3). Also, $y_{{\mathrm{P}}_2,\mathrm{IC}}=1$, $y_{{\mathrm{P}}_2,\mathrm{OC}}=1/2$ (see Section 2.4), and $Y_{\mathrm{IC},\mathrm{OC}}={\mathrm{Y}}_{\mathrm{OC},{\mathrm{O}}_2}=2$ (see Section 2.3.3). Altogether

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{P}}_1& =& (\eta g_1(\mathrm{IC},{\mathrm{O}}_2;\nu )-D){\mathrm{P}}_1,\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{P}}_2& =& (g_2(\mathrm{OC},{\mathrm{O}}_2)-D){\mathrm{P}}_2,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{IC}& =& -g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1+g_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2+D({\mathrm{IC}}_0-\mathrm{IC}),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{OC}& =& (1-\eta )g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1-2g_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2-D\mathrm{OC},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{O}}_2& =& \left(\frac{3}{4}+\frac{17}{40}\eta \right)g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1-\frac{1.9}{4}{g}_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2+D({\mathrm{O}}_{2,0}-{\mathrm{O}}_2).\hfill \end{array}$$

(25)

We track two key quantities, carbon and electrons, through the system. Set

$$\mathrm{C}={\mathrm{P}}_1+{\mathrm{P}}_2+\mathrm{IC}+\mathrm{OC}$$

to be total Cmole concentration in the chemostat to obtain

$$\frac{d}{dt}C=D({\mathrm{IC}}_0-\mathrm{C}),$$

(26)

with solution $\mathrm{C}\left(t\right)={\mathrm{IC}}_0+\mathrm{C}\left(0\right)e^{-Dt}$. So, after a chemostat turnover time $D^{-1}$ or so, $C\left(t\right)$ approaches the constant value $\mathrm{C}={\mathrm{IC}}_0$, the inflow Cmole concentration, to exponentially small error in time. Effectively, thus, the chemostat conserves total Cmoles. Similarly, set the total degree of reduction (DoR) of the system to be

$$\begin{array}{ccc}\hfill \mathrm{DoR}& =& (4.7\mathrm{emole}/\mathrm{Cmole}){\mathrm{P}}_1+(4.1\mathrm{emole}/\mathrm{Cmole}){\mathrm{P}}_2\hfill \\ & & +(0.0\mathrm{emole}/\mathrm{Cmole})\mathrm{IC}+(3.0\mathrm{emole}/\mathrm{Cmole})\mathrm{OC}-(4.0\mathrm{emole}/\mathrm{mole}){\mathrm{O}}_2,\hfill \end{array}$$

Note that

$$\frac{d}{dt}\mathrm{DoR}=-D(4{\mathrm{O}}_{2,0}+\mathrm{DoR}),$$

(27)

with solution $\mathrm{DoR}\left(t\right)=-4{\mathrm{O}}_{2,0}+\mathrm{DoR}\left(0\right)e^{-Dt}$, and thus, after a chemostat turnover time or so, $\mathrm{DoR}\left(t\right)=-4{\mathrm{O}}_{2,0}$, the degree of reduction of the inflow, to exponentially small error in time. Hence, effectively, the chemostat conserves DoR.

In a chemostat, degree of reduction (at least as calculated here) is dictated by the inflow environment, since all reactions conserve it in detailed balance and since all material flows out of the chemostat at the same rate. This contrasts with a biofilm or sparged system where insoluble, soluble, and volatile material may leave the system at different rates. Note that a biofilm can thus have some local control over DoR.

3. Results

3.1. Single Species Chemostat Community

To begin, we consider first the case of a chemostat community of phototrophs only, i.e., ${\mathrm{P}}_2\left(t\right)={\mathrm{P}}_2\left(0\right)=0$. Note that the complementary case of a community of heterotrophs only, i.e., ${\mathrm{P}}_1\left(t\right)={\mathrm{P}}_1\left(0\right)=0$, is not sustainable: ${\mathrm{P}}_1=0$ has the consequence that soluble organic carbon OC is not produced which results in $\mathrm{OC}\left(t\right)\to 0$ which, in turn, results in $g_2\left(t\right)\to 0$ and hence, from (2), ${\mathrm{P}}_2\left(t\right)\to 0$.

In the phototrophic (only) community case, Equations (21)–(25) reduce to

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{IC}& =& -g_1{\mathrm{P}}_1+D({\mathrm{IC}}_0-\mathrm{IC}),\hfill \end{array}$$

(28)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{OC}& =& (1-\eta )g_1{\mathrm{P}}_1-D\mathrm{OC},\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{O}}_2& =& f\left(\eta \right)g_1{\mathrm{P}}_1+D({\mathrm{O}}_{2,0}-{\mathrm{O}}_2)\hfill \end{array}$$

(30)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{P}}_1& =& (\eta g_1-D){\mathrm{P}}_1,\hfill \end{array}$$

(31)

with $f\left(\eta \right)=3/4+(17/40)\eta$ being the net yield of oxygen per carbon fixed, see (10). The coefficients 3/4 and 17/40 arise from degree of reduction details. Note that the equation order has been changed from earlier; the population equation is now listed after the chemical concentration equations for reasons of convenience in the following. The first term in (28) measures usage rate of inorganic carbon in photosynthesis, which produces new biomass (first term of (31)) and soluble organic carbon (first term of (29)) as well as oxygen, some of which is consumed, however, in the production of soluble organic carbon (first term of (30)). Terms involving D measure rates of wash in or out of the chemostat. Note that organic carbon ($\mathrm{OC}$) decouples from the other quantities—the dynamics of IC, O${}_2$, and P${}_1$ are all independent of OC. Hence, system (28)–(31) is effectively three dimensional. We keep OC, though, because of its importance in the two species community dynamics to follow, and also because of its place in conservation of carbon and of degree of reduction.

3.1.1. Steady States

Our interest is in the role of photorespiration in long time community behavior. As is often the case in chemostat models, long time behavior reduces here to the study of steady state solutions. Equations (28)–(31) have two possible types of steady states: (1) the washout solution ${\mathrm{P}}_1\left(t\right)=0$ with $\mathrm{IC}\left(t\right)={\mathrm{IC}}_0$, $\mathrm{OC}\left(t\right)=0$, and ${\mathrm{O}}_2\left(t\right)={\mathrm{O}}_{2,0}$, which exists for all parameter choices though is not always stable, and (2) the viable solution ${\mathrm{P}}_1\left(t\right)={\mathrm{P}}_1^{\ast }>0$ with $\mathrm{IC}\left(t\right)={\mathrm{IC}}^{\ast }$, $\mathrm{OC}\left(t\right)={\mathrm{OC}}^{\ast }$, and ${\mathrm{O}}_2\left(t\right)={\mathrm{O}}_2^{\ast }$. For a viable solution, (31) requires that

$$\eta \left(\mathrm{IC},{\mathrm{O}}_2\right)g_1\left(\mathrm{IC},{\mathrm{O}}_2\right)=D$$

(32)

have a nonnegative solution $\left({\mathrm{IC}}^{\ast },{\mathrm{O}}_2^{\ast }\right)$ indicating that biomass production rate balance with washout. Also, by combining Equations (28) and (30), a second equation,

$$f\left(\eta \right)({\mathrm{IC}}_0-\mathrm{IC})={\mathrm{O}}_2-{\mathrm{O}}_{2,0}$$

(33)

is obtained relating IC${}^{\ast }$ and O${}_2^{\ast }$. Equation (33), an equilibrium relationship constraining the ratio of surplus O${}_2$ to IC deficit, is a consequence of carbon fixation stoichiometry combined with degree of reduction balance. If (32) and (33) can be solved with non-negative values IC${}^{\ast }$ and O${}_2^{\ast }$ then the remainder of a viable state steady state is given by

$$\begin{array}{ccc}\hfill {\mathrm{P}}_1^{\ast }& =& \eta ({\mathrm{IC}}_0-{\mathrm{IC}}^{\ast }),\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill {\mathrm{OC}}^{\ast }& =& (1-\eta )({\mathrm{IC}}_0-{\mathrm{IC}}^{\ast }),\hfill \end{array}$$

(35)

using (28) and (29).

Thus, existence and uniqueness of viable solutions reduces to existence and uniqueness of solutions to (32) and (33). This appears to provide two conditions for viability; in fact, though, (32) and (33) can be solved under the single condition that $\eta ({\mathrm{IC}}_0,{\mathrm{O}}_2)g_1({\mathrm{IC}}_0,{\mathrm{O}}_2)=D$ has a solution with ${\mathrm{O}}_2>{\mathrm{O}}_{2,0}$, that is, under the condition that the organism-free (P${}_1=0$, $\mathrm{IC}={\mathrm{IC}}_0$) chemostat is capable of supporting growth under its given dilution rate D and ambient oxygen level O${}_{2.0}$. For the particular choices of $g_1$, f, and η made here, either one or no viable solutions exist, depending on choice of environmental conditions IC${}_0$, O${}_{2,0}$, and D. See Appendix C for details.

3.1.2. Stability of Steady States

To characterize stability, we add a small pertubation to a steady state solution and then watch ensuing dynamics. We summarize results here; see Appendix D for details. Generally, any component of a perturbation to a steady state that introduces excess or deficient total carbon or degree of reduction is washed out of the system on the chemostat turnover time scale $D^{-1}$. Thus understanding perturbation dynamics of the four dimensional system (28)–(31) reduces to understanding dynamics on a two dimensional subsystem, in fact a system that can be interpreted as the phototroph flux mode space and is spanned by the vectors

$${\mathbf{EFM}}_1=\left(\begin{array}{c}\hfill -1\\ \hfill 0\\ \hfill 4.7/4\\ \hfill 1\end{array}\right),{\mathbf{EFM}}_2=\left(\begin{array}{c}\hfill -1\\ \hfill 1\\ \hfill 3/4\\ \hfill 0\end{array}\right),$$

that encode the two phototroph elementary flux modes. Recall Figure 3: perturbation of the viable steady state by increasing or decreasing flux through the photosynthesis-driven biosynthesis mode corresponds to perturbation of the viable steady state solution in the direction ${\mathbf{EFM}}_1$ (one Cmole biomass and 4.7/4 Omoles produced per Cmole inorganic carbon consumed) and, likewise, perturbation of the viable steady state by increasing or decreasing flux through the photorespiration mode corresponds to perturbation of the viable steady state solution in the direction ${\mathbf{EFM}}_2$ (one Cmole organic carbon and 3/4 Omoles produced per Cmole inorganic carbon consumed). Stability in this flux mode space will be discussed for particular steady states below.

Washout State (One Species System). The washout state (P${}_1=0$) is stable or unstable depending on sign of the quantity $\lambda =\eta ({\mathrm{IC}}_0,{\mathrm{O}}_{2,0})g_1({\mathrm{IC}}_0,{\mathrm{O}}_{2,0})-D$. If negative then the steady state is stable, i.e., phototrophs cannot invade, while is positive, then invasion can occur. Note that λ is the net intrinsic biomass production rate at inflow conditions. When a small quantity of phototrophs are added to the system, in the unstable case $\lambda >0$ dynamics of the linearized system effectively reduce to exponential growth on the one dimensional space $\eta {\mathbf{EFM}}_1+(1-\eta ){\mathbf{EFM}}_2$, indicating that the linearized growth dynamics occurs, as to be expected, as a combination of the photosynthesis mode and the photorespiration mode weighted by the branching parameter $\eta ({\mathrm{IC}}_0,{\mathrm{O}}_{2,0})$.

Viable State (One Species System). For the viable state (${\mathrm{P}}_1>0$), dynamics are again characterized by the basis formed by the two mode vectors EFM${}_1$ and EFM${}_2$ and in this case are always stable (i.e., perturbations decay) under the assumptions that derivatives with respect to IC are non-negative, derivatives with respect to O${}_2$ are non-positive, and condition (17) holds. That is, the viable state is stable under the conditions that we consider biologically reasonable.

3.1.3. Viability and Light-Limited Ranges

We suppose that the RuBisCO inverse specificity factor ${\gamma }_1$, see Section 2.3.3, is subject to some influence by the organism itself, at least adaptively if not through direct regulation, leading to some control over the branching function η. Recall

$$\eta =\frac{1}{1+{\gamma }_1\frac{{\mathrm{O}}_2}{\mathrm{IC}}}$$

(36)

and note that ${\gamma }_1=0$ would correspond to the extreme of $\eta =1$ (all fixed carbon goes to biosynthesis) and that ${\gamma }_1=\infty$ would correspond to $\eta =0$ (all fixed carbon goes to photorespiration). Increasing inverse specificity ${\gamma }_1$ corresponds to increasing, relatively speaking, RuBisCO oxygen affinity and hence increasing photorespiration rates. So then which factors might determine, or at least influence, the value of ${\gamma }_1$?

We show in Appendix E that, for the single species solution as described above (including the assumption (17)), the choice ${\gamma }_1=0$ is favored in the following sense: for any fixed, positive value of ${\gamma }_1$ and the resulting steady state population P${}_1^{\ast }\left({\gamma }_1\right)>0$, it is in fact the case that $(d/d{\gamma }_1){\mathrm{P}}_1^{\ast }<0$. That is, the autotroph population decreases with increasing inverse specificity factor, see, e.g., Figure 4, left panel, for example. Hence, as a larger affinity factor corresponds to increased photorespiration, in the single species, static chemostat environment the autotrophs are always disadvantaged by photorespiration in terms of total biomass.

However, maximizing biomass is not necessarily the only consideration. Another important factor might be viability range—solar light intensity varies significantly over the course of a day (or a year) so that capacity to efficiently function over a wide range of photon flux intensities may also be valuable. High light can cause damage and hence require extra resources, and thus is desirable to avoid or mitigate. In this context, non-zero inverse specificity has competing impacts. First, larger inverse specificity increases, per unit inorganic carbon, the usage of photosynthetically generated electrons and oxygen, thus decreasing rate of damage. Second, larger inverse specificity diverts more fixed carbon from biomass, thus decreasing growth rate. Note though that decreased growth rate leads to reduced population biomass and hence increased available inorganic carbon—a smaller population can be a healthier one. Altogether, then, photorespiration can be expected to shift upwards in both the lower and upper photon intensity viability bounds. To understand how, see Figure 4 right panel, a central result of this study, which presents results of a number of solutions of the steady state Equations (32)–(35). Computations used parameters as described in Appendix B and in the caption.

Minimum Photon Flux. The lower-most curve in Figure 4, right panel, shows, as a function of ${\gamma }_1$, the minimum photon flux necessary for a viable population. This curve was computed analytically by using condition (32) to determine photon flux ν as a function of ${\gamma }_1$ at ambient inorganic carbon and oxygen levels IC = IC${}_0$ and O${}_2$ = O${}_{2,0}$, the limiting viability concentrations. (It was also checked against a numerical computation of minimum ν for viability as a function of ${\gamma }_1$). Its form is easily understood in terms of the non-dimensional number ${\gamma }_1{\mathrm{O}}_{2,0}/{\mathrm{IC}}_0$ which measures the ratio of likelihoods of O${}_2$ versus IC RuBisCO binding in relation to ambient or near-ambient conditions, The ambient ratio O${}_{2,0}$/IC${}_0$ we use is 0.05 Omole/Cmole, i.e., 20 times more inorganic carbon than oxygen as measured in carbon and oxygen moles. Thus, for ${\gamma }_1$ less than approximately 20 Cmole/Omole, binding site competition is unimportant at ambient conditions and hence no penalty, at least with respect to minimum photon flux for viability, is paid. However, as ${\gamma }_1$ increases beyond 20, the minimum photon flux rapidly increases, see Figure 4, right panel.

To summarize, the bottom solid curve in Figure 4, right panel, is important in that it shows minimum photon intensity for community viability as a function of ${\gamma }_1$. This curve is, roughly, described by two parameters: (1) the photon intensity at ${\gamma }_1=0$, which is determined by details of photosynthesis rate function $g_1$ as well as choice of chemostat turnover rate D, and more importantly (2) the value of ${\gamma }_1={\mathrm{IC}}_0/{\mathrm{O}}_{2,0}$ (right-most dotted vertical line) above which significant increase in photon intensity is required for viability.

A similar discussion applies for the upper-most curve in Figure 4, right panel, which shows as a function of ${\gamma }_1$ the maximum photon flux allowable for a viable population. Again, IC ≅ IC${}_0$, O${}_2$ ≅ O${}_{2,0}$, at the viability boundary so that the viability photon flux upper bound is only weakly dependent on ${\gamma }_1$ for ${\gamma }_1{\mathrm{O}}_{2,0}/{\mathrm{IC}}_0$ noticeably less than 1, i.e., ${\gamma }_1$ noticeably less than about 20 Cmole/Omole.

Light-Limited to Carbon-Limited Transition. The dashed curve in Figure 4, right panel, computed numerically, measures as a function of inverse specificity the boundary light intensity between light-limited (region below the curve) and carbon-limited (region above the curve) intensities. In the carbon-limited region, i.e., where photons are sufficiently abundant so that photosynthesis is limited by access to inorganic carbon rather than light, excess electrons are present leading to photoinhibition (recall (14)). The cross-over from light limitation to carbon limitation occurs when $ϵ=\delta$, i.e., when electron production rate as measured in capability to process inorganic carbon is equal to cellular inflow rate for inorganic carbon. The right asymptote (${\gamma }_1\to max\left({\gamma }_1\right)$), shown as the horizontal line in Figure 4 is well approximated by setting ${\gamma }_1=\infty$ and solving $ϵ=\delta$ at IC = IC${}_0$, O${}_2$ = O${}_{2,0}$, the viability values. More importantly,, though, we can understand the small ${\gamma }_1$ behavior of this curve as follows. If ${\gamma }_1=0$ then $\eta =1$ so $g_1=D$ and thus $D=g_1=ϵ=\delta$ which, upon solving for IC and O${}_2$, results in values IC${}_1$, O${}_{2,1}$, the ${\gamma }_1=0$ boundary concentrations. When ${\gamma }_1{\mathrm{O}}_{2,1}/{\mathrm{IC}}_1$ is significantly less than one, i.e., ${\gamma }_1$ noticeably less than ${\mathrm{IC}}_1/{\mathrm{O}}_{2,1}$ photorespiration remains relatively insignificant at cross-over and hence cross-over is only weakly dependent on ${\gamma }_1$. For the parameters used here, ${\mathrm{IC}}_1/{\mathrm{O}}_{2,1}\cong 8\times {10}^{-4}$, see left-most vertical dotted line in Figure 4, right panel. For larger values of ${\gamma }_1$, the range of light-limiting photon intensities expands significantly.

Summary. We summarize Figure 4, right panel, as follows.

• Setting ${\gamma }_1=0$, i.e., turning photorespiration off entirely, results in only a single light intensity with a viable, non carbon-limited steady state population. However, at ambient O${}_2$ and IC concentration levels, competition for RuBisCO binding is insignificant for inverse specificities ${\gamma }_1<{\mathrm{IC}}_0/{\mathrm{O}}_{2,0}$. Hence, from the point of view of population viability at least, there is no penalty for allowing RuBisCO oxygenase activity over this inverse specificity range.
• On the other hand, inverse specificities such that ${\mathrm{IC}}_1/{\mathrm{O}}_{2,1}<{\gamma }_1$ result in significantly enlarged light-limited intensity range, so that large enough inverse specificities may have some advantage.
• Assembled, the inverse specificity interval

  $$\frac{{\mathrm{IC}}_1}{{\mathrm{O}}_{2,1}}<{\gamma }_1<\frac{{\mathrm{IC}}_0}{{\mathrm{O}}_{2,0}},$$

  for parameters used here (based on best approximations in comparison to known data) agrees well with measured values of inverse specificity [11,26,27].

It should be noted that while the upper bound ${\mathrm{IC}}_0/{\mathrm{O}}_{2,0}$ is a function of ambient IC and O${}_2$ levels and is thus is somewhat context-independent at least in the absence of other organisms, the lower bound ${\mathrm{IC}}_1/{\mathrm{O}}_{2,1}$ does depend on specifics of the system like dilution rate D and hence may vary under different conditions. More particularly, IC${}_1$ is found by equating $\delta =D$. For δ as defined here, the resulting concentration IC${}_1$ is given by

$${\mathrm{IC}}_1=\frac{D}{r_{\mathrm{IC}}-D}{K}_{\mathrm{IC}}.$$

(37)

Generally speaking, though, the solution of $\delta =D$ will result in a value of IC${}_1$ as a function of some external rate of transport of IC in comparison to internal, cellular transport mechanisms. Given IC${}_1$ then O${}_{2,1}$ is determined stoichiometrically from (33). Hence the ratio ${\mathrm{IC}}_1/{\mathrm{O}}_{2,1}$ is essentially determined by properties of transport of IC to RuBisCO (relative to the rate of transport of IC into the system), with increased rates corresponding to smaller ratio and hence larger favorable inverse specificity range. Carbon concentration mechanisms, though not included here, might have a similar impact.

Altogether, then, the model suggests that there is possible advantage in the form of redox and oxygen stress control by allowing photorespiration with inverse specificity within the range (37), the lower bound of which is under some internal control. In particular, an expanded range of light-limiting photon intensities may result. This may be important as, typically, photon flux varies considerably over time. (It should be noted that our observation is based on steady state results in a time-independent model, though it seems possible that the idea extends to periodically varying systems). Note that increased photorespiration results in reduced biomass, which may be considered a disadvantage. However, it is in part because of reduced biomass that the range of light-limiting photon intensities increases, as reduction in biomass is accompanied by increase in IC availability.

From the point of view of flux mode modeling, photorespiration provides a sort of rate synchronization mechanism; biosynthesis (left mode in Figure 3) is required to produce biomass, i.e., P${}_1$, at rate dictated by chemostat dilution while photon input is independently, and likely conflictingly, determined by photon inflow rate, both of which are not controlled by the phototrophs themselves (IC input rate can be controlled by the organisms including through varying total biomass). Biomass production rate must match chemostat dilution rate, however, so if the dilution and photon inflow rates differ then, in the absence of a photorespiration mode, excess electrons will be produced leading to damage. Presence of a photorespiration mode (center mode in Figure 3), however, allows some of those excess electrons to be shunted away in the form of reduced OC.

Note that biology-related parameters vary in value between different cyanobacterial species and even within the same species under different environmental conditions, e.g., [29]. We do not see this variability as a critical problem here, however, as our aim is to explore qualitative behavior of community interactions as inverse specificity varies from low (low photorespiration levels) to high (high photorespiration levels), regardless of parameter choices. The forms of the curves in Figure 4 are expected to hold under reasonable choices. In particular we are least confident about choices related excess electron damage, effecting mostly height of the top solid curve, Figure 4 right, and photon usage efficiency, effecting mostly height of the horizontal dotted curve Figure 4 right. The left vertical dotted curve in Figure 4 right, which indicates the approximate value of ${\gamma }_1$ above which photorespiration is significant, is dependent on properties of inflow of inorganic carbon about which we are also relatively uncertain, but because of the log scale used is unlikely, in our view, to move a lot under reasonable choice of parameters.

3.2. Two Species Community

Having explored the effects of photorespiration on steady state phototroph behavior in the one species model, we now add a second species, a heterotroph, in order to see if its addition, despite the resulting (indirect) competition for carbon, can in fact lead to an increase in phototroph steady state biomass. Heterotrophs offer two apparent direct benefits to phototrophs: (1) they use oxygen, thus reducing oxidative stress, and (2) they produce carbon dioxide, thus increasing the local inorganic carbon pool. The price paid is that the cyanobacteria must feed these heterotrophs as they cannot utilize inorganic carbon as a food source. Photorespiration provides a means to do so through production and secretion of soluble organic carbon, thus perhaps providing an additional advantage to its existence. Further, though secretion of organic carbon comes at the price of reduced production of new cyanobacterial biomass, doing so via photorespiration also provides additional control of redox balance through lowering net degree of reduction of the fixed carbon. In this section, then, we consider these combined effects, focusing on steady state cyanobacterial biomass as a metric.

The equations for the two species community are as in (21)–(25), rewritten as

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{IC}& =& -g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1+g_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2+D({\mathrm{IC}}_0-\mathrm{IC}),\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill \frac{d}{dt}\mathrm{OC}& =& (1-\eta )g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1-2g_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2-D\mathrm{OC},\hfill \end{array}$$

(39)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{O}}_2& =& \left(\frac{3}{4}+\frac{17}{40}\eta \right)g_1(\mathrm{IC},{\mathrm{O}}_2;\nu ){\mathrm{P}}_1-\frac{1.9}{4}{g}_2(\mathrm{OC},{\mathrm{O}}_2){\mathrm{P}}_2+D({\mathrm{O}}_{2,0}-{\mathrm{O}}_2),\hfill \end{array}$$

(40)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{P}}_1& =& (\eta g_1(\mathrm{IC},{\mathrm{O}}_2;\nu )-D){\mathrm{P}}_1,\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill \frac{d}{dt}{\mathrm{P}}_2& =& (g_2(\mathrm{OC},{\mathrm{O}}_2)-D){\mathrm{P}}_2.\hfill \end{array}$$

(42)

These equations are the same as the single species ones (28)–(31) except with the addition of source/sink terms proportional to $g_2{\mathrm{P}}_2$ in each of (38)–(40) as well as the new Equation (42) describing heterotroph biomass. In effect we are adding the third elementary flux mode, recall Figure 3, into the system, with all of its component reactions occuring at rate $g_2$.

3.2.1. Steady States and Stability

Equations (38)–(42) have three types of steady state solutions: the washout solution with ${\mathrm{P}}_1={\mathrm{P}}_2=0$, the single species solution with ${\mathrm{P}}_1>0$, ${\mathrm{P}}_2=0$, and the coexistence solution ${\mathrm{P}}_1,{\mathrm{P}}_2>0$. Note that a fourth type of steady state with ${\mathrm{P}}_1=0$, ${\mathrm{P}}_2>0$ is not possible; if ${\mathrm{P}}_1=0$ then the heterotrophs will be washed out of the system. We consider first the washout and single species states, each with P${}_2=0$, and report results of stability analysis here, again referring to Appendix D for details. The coexistence steady state, with P${}_2>0$, is explored numerically later. Note that if P${}_2=0$ then (38)–(41) reduce, essentially, to (28)–(31), so that steady states for the washout and single species systems are the same as previously (with the addition that P${}_2=0$).

Generally, as before, any component of a perturbation to a steady state that introduces excess or deficient total carbon or degree of reduction is washed out of the system on the chemostat turnover time scale $D^{-1}$. Thus, understanding perturbation dynamics of the five dimensional system (38)–(41) reduces to understanding dynamics on a three dimensional subsystem, now spanned by all three of the elementary flux modes, given in vector form by

$${\mathbf{EFM}}_1=\left(\begin{array}{c}\hfill -1\\ \hfill 0\\ \hfill 4.7/4\\ \hfill 1\\ \hfill 0\end{array}\right),{\mathbf{EFM}}_2=\left(\begin{array}{c}\hfill -1\\ \hfill 1\\ \hfill 3/4\\ \hfill 0\\ \hfill 0\end{array}\right),{\mathbf{EFM}}_3=\left(\begin{array}{c}\hfill 1\\ \hfill -2\\ \hfill -1.9/4\\ \hfill 0\\ \hfill 1\end{array}\right),$$

see Figure 3. Perturbation by increasing or decreasing flux through the photosynthesis-driven biosynthesis mode corresponds to perturbation in the direction ${\mathbf{EFM}}_1$ (one Cmole biomass and 4.7/4 Omoles produced per Cmole inorganic carbon consumed) and perturbation by increasing or decreasing flux through the photorespiration mode corresponds to perturbation in the direction ${\mathbf{EFM}}_2$ (one Cmole organic carbon and 3/4 Omoles produced per Cmole inorganic carbon consumed). The new vector EFM${}_3$ corresponds to perturbation that increases or decreases flux through the heterotroph biosynthesis mode (one Cmole biomass and 1 Cmole inorganic carbon produced per two Cmoles organic carbon and 1.9/4 Omoles consumed).

Washout State (Two Species System). The two species washout state (${\mathrm{P}}_1={\mathrm{P}}_2=0$) is, as in the one species washout case, unstable if $\lambda =\eta g_1-D$ is positive and stable if λ is negative. As before, λ is the net intrinsic phototroph biomass production rate at inflow conditions. Also as before, when a small quantity of phototorphs are added to the system, in the unstable case $\lambda >0$, dynamics of the linearized system effectively reduce to exponential growth on the one dimensional space $\eta {\mathbf{EFM}}_1+(1-\eta ){\mathbf{EFM}}_2$, indicating that the linearized growth dynamics occurs, as to be expected, as a combination of the photosynthesis mode and the photorespiration mode weighted by the branching parameter $\eta ({\mathrm{IC}}_0,{\mathrm{O}}_{2,0})$. Note that the heterotroph cannot invade as it requires an already established population of phototrophs (with corresponding finite supply of organic carbon) before it can become viable.

Single Species State: Invasion (Two Species System). We consider for several purposes the single species state (${\mathrm{P}}_1>0$, ${\mathrm{P}}_2=0$). Note that this solution is identical to that in the single species case as discussed in Section 3.1.2 and Appendix C except with the additional component ${\mathrm{P}}_2^{\ast }=0$. We assume that this state is linearly stable to perturbations that do not introduce heterotrophs and ask what happens if a small amount of heterotrophs are added. In other words, how does the otherwise stable heterotroph-free system respond to a perturbation including heterotrophs? This is the invasion problem. In the case that invasion occurs, obviously there is some benefit from the phototrophic population to the invading, heterotrophs as they cannot survive in the chemostat by themselves. Linear analysis can provide some information on the specifics of this advantage.

A key observation here is that the five dimensional system (38)–(42) essentially reduces to the single species, four dimensional one (28)–(31) when P${}_2=0$. In this four dimensional reduced system, dynamics of the single species state are stable, so that the full invasion dynamics are effectively restricted to a complementary one dimensional space. This space is necessarily a linear combination of the three mode vectors EFM${}_1$, EFM${}_2$, and EFM${}_3$. Notably, in the case of large $g_2({\mathrm{OC}}^{\ast },{\mathrm{O}}_2^{\ast })$, instability dynamics are dominated by the heterotrophic growth mode ${\mathbf{EFM}}_3$. When growth is not as dominant, the role of phototroph flux modes in maintaining carbon and DoR balance is more evident. The governing quantity is $\lambda =g_2-D$; if $\lambda >0$ then heterotroph invasion occurs, and if $\lambda <0$ then heterotrophs are unable to invade.

An interesting question here is whether perturbations that include introduction of heterotrophs, i.e., positive perturbation of P${}_2$, result in both successful invasion of heterotrophs as well as, simultaneously, increase in phototroph biomass. We consider this question in the case of large $g_2({\mathrm{OC}}^{\ast },{\mathrm{O}}_2^{\ast })$, see Appendix F for details. Note that, as dynamics are dominantly in the direction of ${\mathbf{EFM}}_3$, then the P${}_1$ component of the perturbation dynamics is small. It is, however, positive as in this case the intuition that addition of heterotrophs, at least initially, increases inorganic carbon concentration and decreases oxygen concentration is correct. Hence, the immediate effects of heterotroph invasion on the phototrophs are mildly positive. Of course, the more important question is of long time effects, which will be considered next.

3.2.2. Two Species Consortium Steady State

We rely on numerical computations to investigate two species steady states. Parameters are as used previously for single species computations with the addition of parameters connected to the heterotrophic biomass mode, see again Appendix B. See Figure 5 for a typical numerical comparison of the single species steady state biomass (as in Figure 4, left panel) and two species steady state biomasses, as functions of inverse specificity. As in Figure 4 left panel, photon intensity is held at a fixed represntative level of 50 $\mathsf{\mu }$E. For the parameters chosen, steady state conditions are carbon-limited for inverse specificity smaller than, approximately, ${10}^{-2}$, and light-limited for larger inverse specificity.

Note that phototroph biomass is larger in the two species community than in the one species community for all values of inverse specificity that allow a viable population. This can be understood as resulting from reconversion of some of the dissolved organic carbon back to inorganic form through heterotroph respiration, recall (19), where it is available for photosynthesis, as opposed to the single species system where all dissolved organic carbon is flushed from the chemostat. Increase in biomass is most noticeable in the interval corresponding to actual measured environmental values of ${\gamma }_1$ where, in the model results, light is limiting and heterotroph biomass is largest

Intuition might suggest that introduction of heterotrophs into a phototoph population would result in increase in IC and decrease in O${}_2$, both as a consequence of respiration. And indeed, such may initially be the case, see the invasion discussion above. However, for later times numerics suggest otherwise near steady state. Dissolved carbon and oxygen, for the same computation, are shown in Figure 6, left and middle respectively. Dissolved inorganic carbon levels are similar for both one and two species communities, with consequently matching protection from photorespiration against high light intensity in both communities. The similarity as well in inorganic carbon concentrations is a consequence of the steady state rate constraint $\eta g_1=D$, see (41); in the small ${\gamma }_1$ carbon limited regime, $\eta \cong 1$ so IC is determined by $\delta =g_1\cong D$ independent of presence or absence of heterotrophs, while in the light limited regime for relatively large ${\gamma }_1$ heterotroph population is low so has little effect and in the light regime with relatively low ${\gamma }_1$, $g_1$ is largely independent of IC and O${}_2$ so that η must be approximately constant, again independently of presence or absence of heterotrophs, see Figure 6 right panel. Note that

$$\left|\frac{{\eta }_{\mathrm{IC}}}{{\eta }_{{\mathrm{O}}_2}}\right|=\left|\frac{{\mathrm{O}}_2}{\mathrm{IC}}\right|.$$

so that in this range, η sensitivity to change in IC is much larger than sensitivity to change in O${}_2$. Hence steady state IC is largely unchanged between the one and two species communities. Dissolved organic carbon, however, is largely absent from the two species community, in contrast to the single species one, as organic carbon is limiting for heterotroph biomass production at all values of inverse specificity and thus is depleted in the two species community. This is perhaps the most dramatic change in chemical environment between the one species and two species environments.

Interestingly, oxygen concentration levels are actually somewhat higher for the two species community, see Figure 6 middle panel, despite oxygen usage via heterotroph respiration. Higher oxygen concentrations can be understood to be a consequence of carbon fixation—increase in reduced carbon in the form of biomass and hence increase in net degree of reduction DoR must be balanced by something. In this model, the only possibility is oxygen. Thus, oxygen concentration increases if biomass does so. This is a generally applicable observation: photosynthetically fixed carbon will be accompanied by more material with low degree of reduction. In an ideal chemostat environment, this balance must be maintained. In other environments, biofilms for example, it may be possible for byproducts like excess oxygen to be transported out of the system while fixed carbon in the form of biomass remains behind.

Note that, in the two species community, as virtually all dissolved carbon goes to biomass for smaller values of ${\gamma }_1$, then in fact two species oxygen levels are approximately independent of inverse specificity ${\gamma }_1$ and this remains so up until the point that ${\gamma }_1$ becomes sufficiently large that significant amounts of inorganic carbon go unutilized. Only then does oxygen concentration significantly decrease as a function of ${\gamma }_1$, though it always remains larger than the corresponding concentration in the one species community.

4. Discussion

The one species model. We observe, in the model, that photorespiration shunts reduced carbon away from biomass production and into dissolved, secreted organic carbon, resulting, at least in a single species oxygenic phototroph population, in three principle effects:

• decrease in population biomass,
• increase in population light tolerance,
• and decrease in oxygen concentration.

The first two are connected through inorganic carbon concentrations: decrease in population results in decrease in inorganic carbon demand resulting in increase in inorganic carbon concentration resulting in reduced inorganic carbon limitation at high light intensities. Decrease in oxygen occurs for two reasons: (1) reduced phototrophic biomass results in reduced oxygen production, and (2) photorespiration product is less reduced than biomaterial, so its production results in less oxygen as a consequence of degree of reduction balance.

The increase in light tolerance and decrease in oxygen concentration suggests an advantage to photorespiration. However, reduction in population size suggests the possibility of fitness deficit in comparison to a population that does not photorespire. We have not modeled such a competition here. However, it should be noted that we impose constant light intensity, and that it is not clear what effects variable light intensity, particularly transient peaks in intensity, might have on a competition of two species, one of which grows more efficiently in low light conditions and the other of which is better protected in high light conditions.

In our set-up, RuBisCO oxygenase activity (which we identify with photorespiration) can serve as a differential of sorts able to synchronize influx of photons with influx of inorganic carbon. Using estimates of parameters, we find an interval of values of the inverse specificity ${\gamma }_1$ which, on the one hand, result in population levels for which, over an increased range of photon intensities, light is limiting but also for which, on the other hand, biomass synthesis is not excessively quenched to the point of reducing the photon intensity range of viability. Though biomass is decreased, the increase in range of “healthy” light intensities might suggest more resilience to light intensity variations, i.e., increased ecological structural stability [30]. The upper bound on ${\gamma }_1$ is related to background inorganic carbon and oxygen concentrations (IC${}_0$ and O${}_{2,0}$ here) and thus may be relatively independent of model details. The lower bound on ${\gamma }_1$ is related to organismal transport rates for inorganic carbon and is perhaps more model dependent, though also allows the possibility of organismal control. In any case, the optimal interval we find for inverse specificities seems to be consistent with measured values over a range of organisms.

The two species model. A steady source of photorespiration-derived organic carbon begs the introduction of a heterotroph population to consume it, so we also modeled a phototroph-heterotroph consortium. Obviously, the heterotroph population benefits from the interaction as it cannot survive in the model chemostat without the organic carbon supplied by the phototrophs. The phototrophs, on the other hand, retain their added tolerance to light intensity but see three principle new effects:

• biomass increase,
• reduction in dissolved organic carbon,
• and oxygen concentration increase.

Biomass increase occurs because of increased inorganic carbon availability as a consequence of heterotroph respiration. Note that organic carbon, here glycolate, could have inhibitory effects, so its consumption by heterotrophs might also potentially increase growth rates, though this effect has not been included in the model. Increase in oxygen concentration is somewhat surprising as heterotrophs consume oxygen during respiration, but occurs again as a consequence of reduction balance: overall increase in reduced biomass must be balanced within the model by increase in oxygen. More directly, the increase in phototroph biomass leads to increased oxygen production. This may be, at least to an extent, an artifact of model simplicity. A more complex model could retain reduction balance through oxidized material other than oxygen. Also, simplicity of the chemostat itself requires that all material be washed out at the same rate whether reduce or oxidized. A more complex system might do otherwise, for example removing dissolved oxygen (e.g., gas sparging that removes O${}_2$) faster than particulate biomass (e.g., biomass fixed in a biofilm).

The phototroph-heterotroph consortium is a more efficient consumer of inorganic carbon than the photorespiring phototroph population alone, and presence of organic carbon suggests that heterotrophs could be expected to join a photorespiring phototroph population. Hence, it may be that the question of competing photorespiring vs. non-photorespiring phototrophs may be the wrong one. Rather, non-photorespiring phototrophs should be asked to compete against a combined photorespiring phototroph-heterotroph consortium.

Connecting flux mode models to population scale models. Mathematically and physically, rates are natural quantities at the flux mode level whereas concentrations, including biomass, are natural at the population and environmental level. We find here that rate functions (in the population model) serve to translate cell scale flux modes into the larger scale population level, where they then determine external concentrations in combination with large scale transport constraints. Flux modes themselves naturally appear, mathematically, in near-steady state dynamics and are relatable to eigenvectors which in turn are natural structures for dynamics. The process of converting flux modes to rate functions is in principle automatable and should be a part of the overall program of extracting information from ‘omics data.

Conversely, the mathematical issues involved in the inverse process of determining how large scale effects influence flux mode regulation are interesting ones and only addressed indirectly here. Generally speaking, microbial communities can have metabolic capabilities available to organisms and to the overall community. This raises the question—how can these capabilities be best deployed to utilize available resources? Rate and stoichiometric constraints still apply, and steady states or, more generally, asymptotic states, can be computed though likely not in a unique way. However, there may be many branching type parameters over which the community has at least some control. Optimality becomes a question of distribution of resource flow (here carbon and electrons) between available pathways in the most efficient manner. Even in the system studied here, with a small number of well defined pathways and a relatively simple physical environment, the effects of that environment on pathway optimization are subtle and influential. The environment imposes rates at steady state and it also determines response to perturbation. These constraints as well as those arising from community ecology may easily be overlooked without considering the physical context of the biological system.