Almost all life on our planet ultimately depends on harvesting the light energy provided by the sun and the subsequent conversion of atmospheric CO2 and other inorganic nutrients into the building blocks of life. As one of the key inventions in evolution, oxygenic photosynthesis has transformed life on Earth and dominates the Earth's primary productivity today (Lane, 2002; Morton, 2009). Beyond their evolutionary and ecological importance, phototrophic organisms are an essential resource for humankind and provide almost all food, feed, and fiber required to sustain human life on this planet with more than 7 billion inhabitants. Many of our strategies to master the challenges of the 21st century will inevitably rely on the growth of phototrophic organisms. Making better use of the sun's light energy while avoiding past mistakes of industrial agriculture related to water usage, energy expenditure, eutrophication, and land use are necessary steps for a sustainable future.

#### *Edited by:*

*Alberto Marin-Sanguino, Technische Universität München, Germany*

#### *Reviewed by:*

*Monika Heiner, Brandenburg University of Technology, Germany Ganesh Sriram, University of Maryland College Park, USA*

> *\*Correspondence: Ralf Steuer<EMAIL_ADDRESS>
#### *Specialty section:*

*This article was submitted to Systems Biology, a section of the journal Frontiers in Bioengineering and Biotechnology*

*Received: 18 February 2016 Accepted: 09 December 2016 Published: 26 December 2016*

#### *Citation:*

*Westermark S and Steuer R (2016) Toward Multiscale Models of Cyanobacterial Growth: A Modular Approach. Front. Bioeng. Biotechnol. 4:95. doi: 10.3389/fbioe.2016.00095*

Phototrophic microorganisms, in particular cyanobacteria, hold great promise as a renewable resource. Cyanobacteria are able to grow with high yield under adverse conditions and their cultivation does not rely on traditional farmland or fresh water. To make use of the biotechnological potential of cyanobacteria, however, requires further understanding of the organization of phototophic growth. While many aspects of phototrophic growth are well known and many details of photosynthetic functioning have been unraveled by decades of active research, it still remains a considerable challenge to understand the individual cellular processes in the context of a living cell.

To this end, the construction of computational models of cellular processes offers the possibility to investigate the emergent properties that arise from interacting processes. Corresponding to the path of experimental research, however, to date almost all computational models involving cyanobacterial functioning and growth focus on the inner workings of individual processes, such as the path of electrons in photosystem II or the functioning of the circadian clock. But phototrophic growth is an organismic property. It is not so much an individual process that gives rise to cellular growth, rather it is the interplay of individual processes that bears reproduction and growth of living cells.

In this contribution, we seek to provide an overview on processes relevant to cyanobacterial growth and summarize the available computational models thereof. Our aim is not encyclopedic, that is, we do not aim for a comprehensive account of all available models. Rather, we seek to focus on representative models that may contribute to our understanding of the functioning of a cell as a whole. Our focus are also not the, albeit important, minutiae of individual processes and models, but rather how they can be collated into a coherent whole. Our starting point is a set of existing computational descriptions of cellular processes and their possible interactions. Our ultimate goal is to describe cellular adaptation, cellular resource allocation, and phototrophic growth in complex environments. Or, as more eloquently put by Neidhardt (1999) already more than 15 years ago: "We must solve the cell. That is, we must do our best to design a computer-based model that can predict overall cell behavior for steady states of growth and for transitions between steady states. The model will at first be crude, inaccurate, and a complete failure at some tasks. With increasing refinement based on additional experimental data, the model should gradually improve. Importantly, the model will guide experimental inquiry by indicating areas of inadequate, insufficient, or incorrect information. Vitally, it is only through such modeling of whole-system behavior—that is, of growth that one will learn how near and how far our knowledge takes us toward understanding the living cell."

Our premise is that sophisticated computational models are already available for many of the processes that underlie phototrophic growth. Modeling their interactions, however, is still no trivial task. First, most models focus on the inner workings of the processes they describe—and therefore often do not describe key variables that govern the interaction with other processes. Second, the various subprocesses and time scales involved in a computational description of phototrophic growth typically require the use of different mathematical and computational concepts, which cannot always be easily reconciled within a single computational description. We seek to summarize these different computational descriptions and aim to highlight common variables and interactions. Importantly, we do not necessarily aim at a single unified model that encompassed all aspects of a growing cell. Rather, we argue for a modular approach—a growing set of models that describe aspects of cyanobacterial growth on different temporal and spatial scales. Depending on the research question, and the temporal and spatial scales involved in this particular research question, different descriptions of cyanobacterial functioning may be chosen—and utilized to derive the emergent properties of cellular growth by putting the parts together.

## 2. MODELING PHOTOTROPHIC GROWTH: AN OVERVIEW

Cellular growth is an organismic process that arises from a coordinated interplay of cellular functions. In the following, we briefly describe the key processes relevant to cyanobacterial functioning and growth in complex environments. An overview is provided in **Figure 1**.

Survival and growth of (most) cyanobacteria begins with the absorption of photons facilitated by large light-harvesting antennae, the phycobilisomes, and chlorophyll *a*. The energy harvested from sunlight drives water splitting at photosystem II (PSII). Electrons, derived from water, are provided to the electron transport chain (ETC) and molecular oxygen is released as the byproduct of photosynthesis.

The ETC consists of a number of large protein complexes, mostly located in the thylakoid membrane. Electrons are transferred along the ETC, ultimately resulting in the regeneration of adenosine triphosphate (ATP) and reduced nicotinamide adenine dinucleotide phosphate (NADPH) as energy carrier and reductant, respectively. The functioning of the photosystems and the ETC are complex biophysical processes and objects of intense research. The respective processes are characterized by fast time scales and transitions between a large number of possible states. While a number of detailed computational models of these processes are available, often with a focus on photosystem II, the respective models typically do not describe regeneration of ATP and NADPH, and hence are not straightforwardly connected to other cellular functions.

The ATP and NADPH regenerated by the photosynthetic light reactions play a crucial role for almost all other cellular processes. Beyond their role as energy donor and reductant, they also serve as important signaling compounds to convey information about the intracellular state. Regenerated ATP and NADPH are utilized to assimilate atmospheric carbon dioxide (CO2). Cyanobacteria possess mechanisms to concentrate inorganic carbon in the vicinity of the CO2-fixing enzyme, the ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO), making use of bacterial microcompartments known as carboxysomes. Compared to the light reactions, the relevant time scales of the so-called dark reactions are significantly slower. Modeling of CO2-concentrating mechanisms (CCMs) typically involves consideration of diffusion and spatial structure.

as well as the circadian clock and its integration into diverse layers of cellular regulation.

The carbon assimilated by the enzyme RuBisCO serves as a substrate to synthesize new cell components, including storage compounds and substrates for cellular respiration. Cellular metabolism involves several hundreds of biochemical reactions, catalyzed by enzymes, as well as spontaneous interconversions, transport, and diffusion processes. From a computational perspective, a description of cellular metabolism must involve several spatial and temporal scales. Detailed models for the action of individual enzymes, such as RuBisCO, exist, based on detailed enzyme mechanisms and elementary reaction steps. Pathways are typically described by combining, often approximative, kinetics of the involved enzymes into larger kinetic models, described by ordinary differential equations (ODEs).

While current kinetic models of metabolism rarely involve more than a few dozen compounds, cellular metabolism is increasingly analyzed using large-scale metabolic reconstructions and constraint-based computational methods. Metabolic network reconstructions are based upon the predicted gene content deduced from genomic DNA and aim to provide an unbiased and comprehensive account of all interconversions of small molecules inside a single cell or a compartment. Unlike kinetic models, metabolic reconstructions only make use of the stoichiometric properties of the respective interconversions, and manually curated reconstructions have been reported for a number of cyanobacteria (Knoop et al., 2010, 2013; Montagud et al., 2010; Nogales et al., 2012; Saha et al., 2012; Vu et al., 2012; Yoshikawa et al., 2015). Highly efficient computational methods exist that allow for the analysis of networks that consist of several hundreds of biochemical reactions and other molecular interconversions. These computational methods, however, are challenging to reconcile and integrate into more traditional enzyme kinetic models of metabolism.

In addition to the photosynthetic light reactions and cellular metabolism, cyanobacterial functioning also involves a large number of regulatory processes. Most prominent is the cyanobacterial circadian clock. Unique among all known prokaryotes, cyanobacteria possess a true circadian clock, a self-sustained oscillator that is entrained to an external *zeitgeber*. Since its discovery in the late 1980s, the cyanobacterial circadian clock has been an object of intense research (Pattanayak and Rust, 2014). Early research, however, was mostly focused on the inner workings of the clock, the molecular details of the core clock, and its input pathways. Only recently, interactions between the clock and metabolism and the question how the clock functions within a broader cellular context have been addressed in more detail (Pattanayak and Rust, 2014; Diamond et al., 2015; Shultzaberger et al., 2015).

Correspondingly, a number of quantitative computational models exist that describe the mechanistic details of the cyanobacterial clock, as well as its entrainment to environmental cues—but as yet only few of these models allow for a straightforward integration into a broader cellular context. While there is increasing evidence how the clock is influenced by, and itself influences, photosynthetic light reactions and metabolism, via sensing metabolic activity (Pattanayak et al., 2015) and redox state (Kim et al., 2012) and controlling transcription regulation, the precise evolutionary role of the clock remains insufficiently understood. Elucidating how the circadian clock interacts with other cellular processes and to integrate models of the cyanobacterial circadian clock into a broader cellular context, with the aim to understand how timing mechanisms affect cellular fitness, is a timely question for further computational research.

Energy metabolism and growth, like all cellular processes, are also dependent on the transcriptional and translational machinery. The transcriptional landscape of commonly cultivated cyanobacteria, such as *Synechocystis* sp. PCC 6803, is increasingly known (Kopf et al., 2014) and a number of studies investigated transcriptional rhythms in the presence of light–dark phases (Lehmann et al., 2013; Beck et al., 2014). Of particular interest is also the role of small regulatory RNA (sRNA) to coordinate cellular processes. For several cyanobacterial strains, most notably *Synechocystis* sp. PCC 6803, substantial sRNA transcription, intragenic transcripts, and antisense transcripts have been reported (Kopf and Hess, 2015). While generic computational models for various possible role of regulatory RNA exist (Legewie et al., 2008), these are currently not integrated within larger computational efforts to understand growth properties and adaptation of cyanobacteria.

As a key driver of cellular functioning and growth, global gene expression is believed to be under direct control of the circadian oscillator, mediated by the topological properties of the cyanobacterial chromosomes. It was shown that the superhelicity of the DNA undergoes rhythmic changes that drive global changes in gene expression (Woelfle et al., 2007; Vijayan et al., 2009). It is further known that the rate of transcription also depends on the local supercoiling status of DNA. Vice versa, supercoiling depends on the cellular energy status, since the extent of supercoiling achieved by the DNA gyrase is strongly dependent on ATP hydrolysis. For heterotrophic organisms, specifically *Escherichia coli*, these observations have led to the proposal of homeostatic control and a feedback loop between the intracellular ATP/ADP ratio, DNA supercoiling, transcription, and again changes in the ATP/ADP ratio (Wijker et al., 1995). Closely related ideas have been put forward in the context of ultradian rhythms in yeast where global partitioning of anabolism and catabolism might be mediated by ATP feedback loop on chromatin architecture (Amariei et al., 2013). In the case of cyanobacteria, a global feedback between cellular energy state, DNA supercoiling, and transcription might mediate between global transcription rhythms, the light reactions as the source of cellular energy, and the circadian clock. Such global feedbacks are currently not explicitly considered in models of cyanobacterial growth and are challenging to implement because of the diverse layers of cellular regulation involved.

Parallel to the efforts of molecular biology to understand the mechanistic and biophysical basis of the processes involved in phototrophic growth, there is a rich history of phenomenological phytoplanktonic growth models. Phenomenological growth formulations are typically employed in models of marine ecosystems and food webs, as well as biogeochemical models to understand the global response of ecosystems to environmental changes. Phenomenological growth models often employ Monod-type equations to describe uptake of a limiting nutrient. Following the early work of Droop (1968), also more sophisticated approaches exist to describe variable internal quotas, see Droop (1983) for an overview. The dynamics of simple phytoplankton growth models are typically based on empirical parameter fitting, rather than an outcome of the underlying cell physiology, and involve strong simplifications, such as using a constant carbon-to-nitrogen (C:N) stoichiometry and absence of photoacclimation (Ayata et al., 2013). It has been pointed out that a major shortcoming of such models is their limited ability to produce true emergence in marine ecosystem models (Allen and Polimene, 2011). Specifically, these models do not evolve to new states not already incorporated in their formulation that makes them unsuitable to properly predict ecosystem changes under changing environmental conditions. As argued by Allen and Polimene (2011), the path forward is to place more emphasis on the underlying intracellular processes—resulting in physiological growth formulations that allow for trade-offs between resource allocations of physiological activities, and hence the possibility to produce biogeochemical and ecological dynamics as emergent properties. Preliminary models, albeit still limited, that combine a detailed description of photosynthesis and phytoplankton growth are already available (Kroon and Thoms, 2006).

In the following, we seek to discuss selected computational models related to cyanobacterial functioning and growth in more detail. Our view is that cyanobacterial physiology depends on interacting cellular processes that can be interpreted as functional "modules", such as the photosynthetic light reactions and the ETC, carbon uptake mechanisms, cellular metabolism, the circadian clock, as well as the transcriptional and translational machinery and its regulation. For many of these modules, reasonable computational descriptions already exist, whereas other processes, for example, the coordination of cell cycle events in relation to metabolism (Asato, 2005, 2006), have not yet been subject to computational studies.

Our aim is to highlight the common variables and known interactions between the processes relevant to cyanobacterial functioning and growth. In this respect, a particular challenge is the wide range of computational approaches and methods used. Models of cellular processes may take many forms, ranging from spatial versus non-spatial, stochastic versus deterministic, population level versus single cell level, and continuous versus discrete descriptions. See **Figure 2** for an overview. Notwithstanding the technical challenges, we believe that the integration of different aspects of cellular growth, and their respective computational representation, is a prerequisite toward understanding the living cell. We seek to understand how phototrophic growth functions and how it is regulated. How does the coordination of physiological functions work in order to synthesize the right macromolecules at the right time? Which

FIGURE 2 | Models of cellular processes are highly diverse and may involve a wide range of computational concepts and methodologies. At the core of the modeling process is a translation of a biological processes into a formal (mathematical) language. Once this translation is established, the model can be interrogated using the tools of mathematical and computational analysis. The most prevalent representations of cellular processes described in this contribution make use of deterministic ordinary differential equations (ODEs) to describe the time-dependent dynamics of continuous intracellular concentrations, typically on the population level. In the following, such models are denoted as kinetic models and may either make use of heuristic approximate rate equations or rate equations derived from explicit biochemical mechanisms. Models of CCMs typically involve a spatial component and originate from a description based on partial differential equations (PDEs). Models of the light reactions frequently describe transitions between discrete states that occur with a certain state-dependent probability. Flux balance models consists of a set of linear relationships (linear inequality constraints) between variables and make use of linear programing (LP), a method to identify the optimum of a linear objective function. For a more detailed overview on model types, see also Steuer and Junker (2009).

level of detail is required to describe cellular growth? What are the variables and time scales involved?

# 3. MODELS OF THE PHOTOSYNTHETIC LIGHT REACTIONS

Phototrophic growth begins with the absorption of light and its conversion into chemical energy. Despite a number of open questions and the need for further research, many of the fundamental properties of oxygenic photosynthesis have been elucidated in the past century. Owing to the fact that cyanobacteria are the evolutionary ancestors of modern-day chloroplasts, the organization of their photosynthetic ETC is essentially identical to that in algae and green plants (Vermaas, 2001).

In most cyanobacteria, light harvesting is facilitated by large antenna complexes, the phycobilisomes. Phycobilisomes are attached to the cytoplasmic surface of the thylakoid membrane (Mullineaux, 2014). The detailed composition of phycobilisomes is strain specific and depends on light quality, denoted as complementary chromatic adaptation. The energy absorbed by the phycobilisomes is transferred to either photosystem II or photosystem I, or dissipated as heat or fluorescence. The protein complexes of the photosynthetic ETC are embedded within the thylakoid membrane. The key proteins complexes responsible for photosynthetic electron transport are Photosystem II (PSII), the Cytochrome b6f complex (Cytb6f), Photosystem I (PSI), and ATP synthase (ATPase). See **Figure 3** for an overview.

PSII splits water and reduces the plastoquinone (PQ) pool. The latter mediates the transport of electrons from PSII to Cytb6f. At Cytb6f, electrons are transferred to a soluble electron carrier on the luminal side of the thylakoid membrane, either plastocyanine (PC) or cytochrome-*c* (cyt-*c*). At PSI, electrons are transferred to ferredoxin and eventually to NAPDH using light-induced excitation of the PSI reaction center (linear electron transport, LET). Alternatively, electrons from the excited PSI state can be transferred back to PQ and Cytb6f (cyclic electron transport, CET), details of CET are still under debate and insufficiently understood. Photosynthetic electron flow results in a protein gradient across the thylakoid membrane that drives regeneration of ATP by the ATPase.

Unique to cyanobacteria, as opposed to plants and microalgae, is the combination of oxygenic photosynthesis and respiration in the same membrane system using intersecting ETCs and common components (Vermaas, 2001). The respiratory ETC involves the succinate dehydrogenase (SDH), the NADPH dehydrogenase (NDH-1), and terminal oxidases. The PQ pool, the Cytb6f complex, and PC (or cyt-*c*) as soluble electron carrier are involved in respiratory as well as photosynthetic electron transport. While photosynthesis exclusively takes place in the thylakoid membrane, a rudimentary respiratory chain is also present in the plasma membrane (Schultze et al., 2009).

From a computational perspective, photosynthesis in cyanobacteria and microalgae can be described on different levels of complexity. Basic models are closely related to overall growth models in ecology—and typically to reproduce the production of oxygen and the photosynthesis–irradiance (PI) curve of cyanobacteria and microalgae. Early models were derived by Crill (1977), Megard et al. (1984), Eilers and Peeters (1988), and Zonneveld (1998) among others. These models make use of a highly simplified photosynthetic factory or photosynthetic unit (PSU) that encompasses PSII, PSI, and the ETCs. See **Figure 4** for an example. The resulting differential equations for the dependency of photosynthesis on light intensity can often be solved analytically, with a solution analogous to the Haldane equation—an enzyme kinetic equation that was derived for substrates with inhibitory effects at high concentrations. Simple three-state models are suitable to describe basic features of photoinhibition and the PSII repair cycle (Tyystjärvi et al., 1994). In later iterations, the parameters of the basic three-state model were augmented with a more mechanistic interpretation (Han, 2001, 2002), and the models were extended to describe the effects of intermittent light (Rubio et al., 2003). Recently, a basic three-state model was also applied to describe the kinetics of non-photochemical quenching (NPQ), induced by an orange cartenoid protein (OCP), in cyanobacteria (Gorbunov et al., 2011). To this day, simple three-state models remain relevant

to describe overall photosynthetic activity, in particular in bulk models to assess productivity in photobioreactors and related industrial application (Nedbal et al., 2010; Bernard, 2011). With respect to understanding interactions between cellular processes, a drawback of highly simplified growth models is their insufficient representation of intracellular parameters, such as no explicit PQ pool, no explicit regeneration of ATP and NADPH, and lack of alternative electron transport.

Beyond overall bulk models of photosynthesis, there is a significant history of biophysical models to understand oxygen evolution and chlorophyll fluorescence transients, often with a focus on PSII, as well as to understand specific properties, such as energy distribution in the photosynthetic apparatus (Butler and Strasser, 1977; Butler, 1978) or, more recently, excitation transfer in the PSII membrane (Amarnath et al., 2016). Early kinetic models were described by Mar and Govindjee (1972), a more elaborate model was put forward by Holzwarth et al. (2006) and later analyzed by Nedbal et al. (2007). Further elaborate models of this kind were developed by Lazár (2003) and Zhu et al. (2005). The former was refined and extended by Jablonský and Lazár (2008); different approaches were later compared by the same authors (Lazár and Jablonský, 2009). Common to these models is a focus on chlorophyll fluorescence emission, and to a lesser extent oxygen evolution, as the main output variables. While relevant for biophysical research, the respective models cannot be straightforwardly integrated into more comprehensive models of phototrophic growth, due to the focus on fast time scales and specific output variables. We note that the interpretation of results obtained from pulse-amplitude modulated (PAM) fluorimetry significantly differs between cyanobacteria and plants (Schuurmans et al., 2015; Acuña et al., 2016), with modeling approaches focusing almost exclusively on the latter.

Models that explicitly describe the photosynthetic electron transport chain and subsequent reactions, in addition to PSII, are more suitable to integrate into the context of a living cell. To this end, a small number of models exist (Berry and Rumberg, 2000; Vershubskii et al., 2014), typically based on ODEs. An elaborate model of this type was proposed by Laisk et al. (2006), developed to understand the photosynthetic process from light absorption to sucrose synthesis. The model neglects many of the detailed biophysical properties of earlier models (Zhu et al., 2005), such as an explicit representation of the s-states that describe the cyclic reactions of the oxygen-evolving complex (Kok et al., 1970). The model instead provides a combination of whole chain electron transport and carbon assimilation processes, including nonphotochemical quenching, chlorophyll fluorescence, and (albeit

FIGURE 4 | The photosynthesis–irradiance (PI) curve for a minimal model of photosynthesis (Han, 2002). The ETC is described by a photosynthetic unit (PSU) that exists in an open or reactive state. After being subjected to light, the PSU transits to a closed or activated state (PSU\*). Excessive absorption results in photodamage and an inhibited state PSU*<sup>d</sup>* with rate constant *kd*. (A) Using simple ODEs based on mass action kinetics results in typical PI curves. The overall functional form is similar to the Haldane equation, an equation derived for enzymes whose substrates have an inhibitory effect at higher concentrations. (B) The reaction scheme. Similar models can describe non-photochemical quenching in which the damaged state corresponds to a quenching state (indirectly) activated by light (Gorbunov et al., 2011). All values are reported in arbitrary units (a.u.). The ODEs used to generate the figure are provided in the Appendix.

simplified) photorespiration. Using a similar approach, Zhu et al. (2013) described a detailed dynamic model of leaf photosynthesis, based on ODEs, from light capture to carbon assimilation, and incorporates the previous partial model of the same authors (Zhu et al., 2007) augmented by additional reactions of the ETC. Both models focus on C3 plant metabolism, but similar approaches are feasible for cyanobacteria.

Importantly, both models provide a sufficient level of detail to interface with other cellular processes and include ATP regeneration and reduction of NADPH, photorespiration, alternative electron transport, as well as an explicit representation of the PQ pool and lumenal pH. Owing to the focus on plant C3 metabolism, neither of the models describe peculiarities of cyanobacteria, such as shared components between the photosynthetic and respiratory ETC and the resulting differences in regulation.

Selected models of the ETC and the photosynthetic light reactions are summarized in **Table 1**. Main challenges for the development of corresponding models for cyanobacteria are to incorporate the respiratory ETC, as well as to incorporate the specific alternative electron sinks of cyanobacteria. As an interface to other cellular processes, discussed below, a model of the cyanobacterial ETC should include regulatory switches in cyanobacterial photosynthesis, regulation of light harvesting including regulation of the orange carotenoid protein, state transitions that control the relative energy transfer from phycobilisomes to PSII versus PSI, and alternative electron sinks that serve as "electron valves" and prevent overreduction of the ETC, among other features that are relevant for cyanobacterial functioning and growth (Mullineaux, 2014). Relevant exchange variables are ATP and NADPH regeneration, the state of the PQ pool, leakage of reactive oxygen species (ROS), and oxidation of metabolites for cellular respiration.

#### TABLE 1 | Selected models of the photosynthetic light reactions and the electron transport chain (ETC) in plants and cyanobacteria.

*This table and all following tables are not exhaustive but highlight selected models that are of particular relevance to describe aspects of cyanobacterial growth. Model details, such as the exact number of variables, are listed only if they can be unambiguously sourced from the original publication.*

# 4. KINETIC MODELS OF CELLULAR METABOLISM

The energy harvested by the photosynthetic light reactions drives the assimilation of inorganic carbon and the synthesis of storage compounds and building blocks for cellular growth. Photoautotrophic metabolism involves the uptake of inorganic carbon facilitated by CO2-concentrating mechanisms (CCMs), assimilation of CO2 by RuBisCO, and the subsequent synthesis of cellular building blocks mediated by a network of metabolic reactions. Computational concepts used to describe cyanobacterial metabolism have been discussed previously (Steuer et al., 2012), here we focus on the integration of such descriptions into integrative models of phototrophic growth. In particular, models of metabolism are highly diverse and span multiple orders of magnitude with respect to the time scales and number of variables involved. As highlighted previously (Steuer and Junker, 2009), no single universal computational methodology exists that is suitable to describe all relevant aspects of metabolic functioning. Rather, a hierarchy of computational approaches exists, ranging from detailed kinetic models of individual enzymatic reactions, based on ODEs, to large-scale stoichiometric reconstructions that are evaluated using constraint-based analysis.

The basic building blocks of metabolism are the actions of individual enzymes and their respective reaction mechanisms. Computational modeling of enzyme kinetics is well understood (Steuer and Junker, 2009; Sauro, 2014), even though specific reaction mechanisms and atom transition maps are not yet comprehensively available and must be confirmed on an individual per-reaction basis. Detailed kinetic models of key reactions have been proposed in the literature, most notably for RuBisCO, the key enzyme of the Calvin–Benson cycle (Witzel et al., 2010). Following the rules for overall rate equations of enzyme kinetic mechanisms, multiple reaction steps can be integrated into larger models of cellular pathways. Corresponding detailed kinetic models of metabolic pathways exist since the late 1950s, see Steuer and Junker (2009) for a review, and several kinetic pathway models have since been proposed for phototrophic plant metabolism. Of particular interest is the Calvin–Benson cycle and the adjacent carbon metabolism. Among the first computational descriptions of phototrophic C3 carbon metabolism were the models proposed by Milstein and Bremermann (1979) and Hahn (1984). The former involves 17 first-order ODEs and 22 parameters. The latter involves 19 state variables and describes the dynamics of Calvin–Benson cycle intermediates, as well as parts of sucrose and starch metabolism. The model was later extended to include photorespiration (Hahn, 1987), and simplified representations were considered (Hahn, 1991). The latter analysis demonstrated that also smaller models are able to reproduce the observed dynamics. Both models, as well as the model of Laisk and Walker (1986), are largely based on mass action kinetics, rather than derived enzyme kinetic equations for kinetic mechanisms. Parallel to the development of detailed kinetic models of core carbon metabolism, a number of biochemical models of photosynthetic CO2 assimilation have been developed that focus on plant-specific properties, such as gas exchange and stomatal conductance, often also incorporating aspects of carbon reduction, see Farquhar et al. (1980)

for an influential early example. Likewise, a significant number of models were developed to investigate the origin of photosynthetic oscillations (Giersch, 1986; Laisk and Walker, 1986; Laisk et al., 1989; Rovers and Giersch, 1995), see Roussel and Igamberdiev (2011) for a recent review.

The prototype of most current enzyme kinetic models of the Calvin–Benson cycle was proposed by Pettersson and Ryde-Pettersson (1988). The model is based on mechanistic non-linear enzyme kinetic rate equations, implemented as ODEs, together with equilibrium mass-action ratios. The model describes the dynamics of the Calvin-Benson cycle under conditions of light and CO2 saturation. The parameterization of the model involved approximately 50 kinetic parameters, sourced from the literature across several plant species. The model of Pettersson and Ryde-Pettersson (1988), like many current and past kinetic models, is therefore not necessarily a model of a single plant species, but must be considered as a prototype model that describes several generic aspects of plant leaf C3 metabolism. The model was later adapted to investigate further aspects of metabolic regulation in phototrophic metabolism (Poolman et al., 2000, 2001), and extended by Zhu et al. (2007) to investigate the reallocation of enzymes of photosynthetic carbon metabolism with respect to optimal nitrogen and protein investment. These models also served as a blueprint for the first detailed kinetic models of cyanobacterial core carbon metabolism. Jablonský et al. (2013) proposed a modified version of the model of Zhu et al. (2007), adapted to describe the cyanobacterium *Synechococcus elongatus* PCC 7942, to investigate the functional consequences of isoenzymes. The majority of parameters were retained from the original models. The model was later refined (Jablonský et al., 2014) to explain the metabolic regulation of primary carbon metabolism, also incorporating transcriptional data as a constraint for model dynamics. Most recently, a kinetic model of the central carbon metabolism of the cyanobacterium *Synechocystis* sp. PCC 6803 was developed to investigate the role of isozymes on metabolic network homeostasis with respect to changes in gene expression induced by different CO2 conditions (Jablonský et al., 2016). In particular, a comparison of model properties indicated that the higher number of isozymes present in the *Synechocystis* sp. PCC 6803 genome compared to the (smaller) genome of *Synechococcus elongatus* PCC 7942 may correspond to a shift of metabolic regulatory strategies from transcriptional control in latter toward post-transcriptional control in the former (Jablonský et al., 2016).

From computational perspective, the kinetic models considered above share several features relevant to the integration into multiscale models of phototrophic growth. In each case, the dynamics of the concentrations of metabolic intermediates are described by ordinary differential equations (ODEs). Rate equations are derived from enzyme kinetic mechanisms, and implemented using (usually reversible) non-linear Michaelis– Menten type functions. The rate equations consider allosteric regulations, as well as other post-translational mechanisms, as far as such interactions are known. The light reactions are typically highly simplified. In the model of Pettersson and Ryde-Pettersson (1988) and its later adaptations, ATP is provided by a single overall reaction (an ATP synthetase) that can be modulated according to light intensity. The concentrations of NADP<sup>+</sup> and NADPH are assumed to be constant. Likewise, enzyme amounts are represented by external parameters, the respective values are part of the maximal reaction velocities Vmax. **Table 2** lists selected kinetic models of central carbon metabolism and the Calvin–Benson cycle. We conjecture that such models provide a reasonable account of metabolite dynamics on short and medium time scales (minutes to few hours) and to metabolic adaptations to brief periods of darkness. As yet, the construction of kinetic models to adequately represent changes in day/night metabolism remains a considerable challenge.

The key factors to integrate models of core carbon metabolism into overall models of phototrophic growth requires an explicit representation of energy (ATP) and redox state (NADPH/NADP<sup>+</sup>) as dynamic variables that allow coupling to the ETC. In this respect, first steps have been made for plant metabolism. The models of Laisk et al. (2006) and Zhu et al. (2013) both integrate the photosynthetic light reactions with a detailed representation of the core C3 carbon metabolism, including photorespiration. Such models provide a framework for guiding engineering efforts and allow for a description of photosynthesis and carbon fixation in response to, for example, changes in photon flux density.

Beyond the integration of light capture and carbon metabolism, significant challenges remain. Allosteric post-translation regulation only covers small to medium time scales. No current kinetic model provides a description of diurnal changes in metabolism and the corresponding switch from carbon assimilation to the mobilization of storage compounds. In addition to redox regulation, such switches are likely to require the inclusion of additional hierarchies of cellular regulation, in particular transcription and possibly regulation by the circadian clock. Switches in metabolic modes are of particular relevance for cyanobacterial metabolism and growth, as they lack the compartmentation of eukaryotic algae and plants. Likewise, current models focus on carbon metabolism and its regulation, limitation of other macronutrients like phosphorus or nitrogen is typically not considered. Nonetheless, in particular with respect to nitrogen, kinetic models can be expected to provide insight into the role of certain metabolites, such as 2-oxoglutarate (2-OG), as signaling compounds (Fokina et al., 2010).

*The models of Laisk et al. (2006) and Zhu et al. (2013) integrate the photosynthetic light reactions and the core carbon metabolism and are already listed in Table 1. Minimal models that investigate oscillatory mechanisms, and models that primarily focus on plant-specific properties, such as stomatal conductance, are not considered.*

A major obstacle for detailed kinetic models also remains the scarcity of enzyme kinetic data. The construction of kinetic models requires detailed knowledge of enzyme kinetic parameters—with typically 4–5 parameters per reaction, including the Michaelis–Menten parameters *K*M for substrates and products, as well as thermodynamic equilibrium values Keq and the specific catalytic activities of enzymes. Even metabolic pathways of moderate size, such as the Calvin–Benson cycle and adjacent reaction, typically consist of 20–30 enzymatic reactions. Therefore, the construction of larger kinetic models, while feasible from a computational point of view, is primarily limited by data availability and data reliability (Srinivasan et al., 2015). To some extent, the scarcity of information about kinetic parameters can be alleviated by explicitly accounting for uncertainty in kinetic models of metabolism—suitable approaches have been proposed recently (Wang et al., 2004; Steuer et al., 2006; Tran et al., 2008; Steuer and Junker, 2009; Murabito et al., 2014) but are not yet widely applied in models of phototrophic growth.

# 5. MODELS OF CARBON-CONCENTRATING MECHANISMS

A characteristic feature of cyanobacterial growth is the use of CO2-concentrating mechanisms (CCMs) to facilitate the uptake and acquisition of inorganic carbon. CCMs allow cyanobacteria to raise the local concentration of CO2 in the vicinity of the carboxylating enzyme RuBisCO, and thereby overcome the comparatively low affinity of RuBisCO for CO2

and depress the competitive oxygenation reaction (photorespiration). Cyanobacterial CCMs typically make use of dedicated microcompartments, the carboxysomes, that separate the assimilation of CO2 by RuBisCO from the rest of the cell. The CCMs of cyanobacteria relies on a number of components. The respective mechanisms are reasonably well understood (Kaplan and Reinhold, 1999; Price et al., 2008; Burnap et al., 2015), see **Figure 5** for a schematic depiction, and the requirement of a quantitative mathematical analysis has recently been highlighted (Mangan and Brenner, 2014). The efficiency of the cyanobacterial CCMs can be characterized by the ratio between the apparent whole-cell affinity for extracellular CO2 and the respective affinity for CO2 of the carboxylating enzyme RuBisCO—with ratios up to 1,000 reported in the literature (Burnap et al., 2015). While many components of the CCMs are constitutively expressed, the expression of specific uptake systems is differentially regulated depending on environmental parameters, in particular light intensity and the availability of inorganic carbon (Ci) (Kaplan and Reinhold, 1999; Burnap et al., 2015).

In addition to their important function to enhance the local CO2 concentration and depressing photorespiration, CCMs may also be involved in dissipating excess light energy (Xu et al., 2008) and might play a role in pH homeostasis. Despite this tight integration with carbon metabolism, none of the current kinetic models of carbon metabolism explicitly accounts for CCMs. Nonetheless, a number of quantitative models of CCMs are available that can be readily incorporated into models of cyanobacterial growth. See **Table 3** for selected computational models of CCMs. Early

transporters and Na+/ HCO3 <sup>−</sup> symporters. The activity of CCMs, including expression of several components, is modulated by environmental parameters, in particular CO2 availability. CO2 leaking from the carboxysomes may diffuse into the medium or is partly converted back to HCO3 <sup>−</sup> (carbon cycling) at the thylakoid membrane using insufficiently understood mechanisms. CCMs utilize cellular energy and might be involved in dissipating excess light energy and play a role in the maintenance of internal pH. Hence, CCMs are integrally tied to cellular metabolism and growth. See Table 3 for selected computational models of CCMs.

#### TABLE 3 | Selected models of cyanobacterial CO2-concentrating mechanisms (CCMs).

*Modeling the CCM typically involves a spatial component. Analysis of the respective partial-differential equations (PDEs) shows, however, that intracompartmental concentration gradients can be neglected. Models of the CCM are therefore typically based on ODEs, rather than PDEs.*

models focus either on simple equations for CO2 and HCO3 − (bicarbonate) flux into and out of the cell (Badger et al., 1985), as well as on arguments based on reaction–diffusion equations (Reinhold et al., 1987, 1991). These models were refined further to include explicit representations of the carboxysomes. Specifically, the model of Fridlyand et al. (1996) considers the various CO2 or HCO3 − fluxes between medium, periplasmic space, cytoplasm and carboxysomes using derived values for geometric parameters, and permeability and diffusion coefficients. The model also considers the energetic consequences of scavenging CO2 that leaks back into the cytoplasm. Models can be adopted to specific cyanobacterial strains, such as the model of Hopkinson et al. (2014) for *Prochlorococcus* spp. MED4. Two recent quantitative models of CCM functioning are available (Clark et al., 2014; Mangan and Brenner, 2014). Both models are based on reaction– diffusion equations that are solved for highly simplified spatial topologies (spherical cells), but otherwise make use of partially divergent assumptions. The focus of Clark et al., 2014 are interspecies differences and a hypothetical carboxysome-free mutant that is of interest in industrial settings with elevated CO2 supply. The model assumes that the carboxysome walls are impermeable to CO2. The model of Mangan and Brenner (2014) assumed that carboxysome permeability is identical for HCO3 − and CO2 and the model explores the range of best parameter values that give rise to a functional and effective CCM. While carboxysome permeability has not yet been measured directly, Mangan and Brenner (2014) concluded that optimal parameter values indeed exist, and transport rates and concentrations derived for these optimal values are in good agreement with known experimental data. Very recently, the model was extended to incorporate the effect of intracellular pH as a key physiological parameter that governs the composition of the Ci pool (Mangan et al., 2016). The "pH-aware" model highlights the utility of quantitative models to evaluate the energetic costs of Ci accumulation for CCMs.

While current models of the CCM do consider the optimality and functioning of CCMs under different intracellular and environmental conditions, they typically do not incorporate explicit models of cellular growth or other cellular mechanisms. However, the high energy demand of Ci transport, either by direct hydrolysis of one ATP per bicarbonate, or indirectly via the costs of ion transport, the costs of synthesizing carboxysome shell proteins, as well as the significant impact of CCMs on the efficiency of carbon assimilation suggest that further integration of models of CCMs into a broader context of cellular functioning is worthwhile to understand the trade-offs and interactions between energy investment, CCM utilization, carbon assimilation, and growth.

## 6. LARGE-SCALE MODELS OF CYANOBACTERIAL METABOLISM

Beyond kinetic models of central carbon metabolism, metabolic networks are increasingly described using large-scale stoichiometric reconstructions (Steuer et al., 2012). Metabolic reconstructions aim to provide a comprehensive account of all possible interconversions of small molecules within a given cell or compartment, including enzymatic reactions as well as non-catalyzed (spontaneous) interconversions, transport reactions, and diffusion. Metabolic reconstructions are derived from annotated genomes with subsequent steps of manual curation and gap filling, see Knoop et al. (2010), Steuer et al. (2012), and Knoop et al. (2013) for recent examples. The description typically involves only knowledge about the stoichiometry of interconversions; knowledge about kinetic parameters (such as Michaelis– Menten parameters) or allosteric regulation is not required.

Nonetheless, large-scale stoichiometric models of bacterial metabolism are highly predictive (McCloskey et al., 2013). The predictive power derives from the fact that the fluxes through enzymatic reactions are not independent. Constraint-based computational methods rely on the fact that under steady-state conditions most intracellular metabolites do not accumulate. The rate of synthesis of any non-accumulating metabolite must therefore approximately equal the rate of consumption of this metabolite. Similar arguments also hold for diurnal metabolism: if, after a full diurnal cycle, the concentration of a given intracellular metabolite is approximately equal to its initial value, then the total flux of synthesis reactions and the total flux of consuming reactions must be approximately equal. The condition of flux balance puts significant constraints on the feasible flux space. Predictions about specific flux solutions are then typically based on assumptions about metabolic optimality. That is, among all feasible flux solutions, constraint-based methods seek to identify a solution that maximizes a given objective function, such as the maximal yield of biomass for a given light input—motivated by the fact that a similar selection might take place during evolution. Predictions from large-scale stoichiometric models are therefore not mechanistic, that is, they are not derived from knowledge about biophysical or biochemical interactions. Rather, predictions are derived from how metabolism *ought* to behave given the assumption that metabolic functioning fulfills certain evolutionary optimality principles.

Computationally, the analysis of large-scale stoichiometric reconstructions is based on methods of linear programing (LP) and is computationally feasible also for networks consisting of several thousands of reactions. The strength of stoichiometric models and constraint-based analysis are questions such as the following: What is the maximal growth yield for a given light or carbon input? Which set of enzymes is essential for the synthesis of certain biomass components? How many distinct biochemical paths exist for the synthesis of certain biomass components and how do these pathways differ with respect to cellular energy expenditure and cofactor utilization? Due to the specific computational methodology, however, a direct integration of large-scale stoichiometric models into kinetic models of metabolism remains challenging. Various extensions toward incorporating dynamics have been proposed (Mahadevan et al., 2002; Kim et al., 2008; Feng et al., 2012; Antoniewicz, 2013), and extensive efforts are undertaken to bridge the gap between kinetic and stoichiometric models (Steuer, 2007; Steuer and Junker, 2009; Chakrabarti et al., 2013; Srinivasan et al., 2015).

 See **Table 4** for selected examples of metabolic reconstructions. Such large-scale reconstructions are valuable tools to derive consistent equations for the (maximal) growth yield with respect to light input, to derive core models of reaction pathways, and to make predictions about maximal product yield in biotechnological applications (Zavřel et al., 2016). In particular, large-scale reconstructions also enable a semiautomated extraction of meaningful core models to facilitate the construction of smaller kinetic models (Erdrich et al., 2015). Analysis of the respective networks, however, is often confined to either a constant light environment, or to heterotrophic growth

*Reconstructions of algae and land plants are not considered.*

on extracellular carbon sources. Only recently stoichiometric analysis of phototrophic metabolism explicitly described different phases of light availability. For example, Knoop et al. (2013) simulated biomass synthesis fluxes over a full diurnal cycle, Muthuraj et al. (2013) used dynamic FBA (dFBA) to capture light-dark metabolism over discretized time intervals, Knies et al. (2015) considered storage metabolites that accumulate and are consumed over a diurnal cycle using a reconstruction of the unicellular alga *Emiliania huxleyi*, Cheung et al. (2014) described a flux balance model that captures interactions between light and dark metabolism in C3 and CAM leaves, and Baroukh et al. (2014) proposed a novel dynamic modeling framework to describe carbon metabolism of unicellular microalgae.

Beyond conventional FBA and related constraint-based methods, there has been increasing interest to evaluate cellular metabolism in terms of a cellular "protein economy" (Molenaar et al., 2009) and to study trade-offs in cellular resource allocation (Goel et al., 2012; Müller et al., 2015)—a theme where, as noted by Schaechter (2015), the bacterial growth physiology of old is connected to the systems biology of today (Stouthamer, 1973; Neidhardt et al., 1990). As one of the first applications involving cyanobacteria, Burnap (2015) formulated a model of autotrophic growth in terms of allocating protein resources among core functional groups, such as the ETC, light-harvesting antennae, and ribosomes. Along similar lines, Rügen et al. (2015) formulated a self-consistent autocatalytic model of phototrophic growth. The model is based on the observation that the macromolecules that constrain cellular growth, including the components of the ETC, metabolic enzymes, and ribosomes, are itself products of metabolism. Phototrophic growth can therefore be formulated as a time-dependent linear optimization problem, such that optimal growth entails a time-dependent allocation of resources during a full diurnal cycle. The approach of Rügen et al. (2015), denoted as conditional FBA, results in dynamic time courses for all involved reaction fluxes, as well as changes in biomass composition over a diurnal cycle. Similar to conventional FBA, models of this kind are not based on mechanistic insight, but rather seek to evaluate the optimality of resource allocation during phototrophic growth. It is expected that methods and applications that go beyond conventional FBA and involve spatial and temporal metabolic modeling based on genome-scale reconstructions of microbial metabolism will play an increasingly important role (Henson, 2015).

### 7. MODELS OF THE CYANOBACTERIAL CLOCK

In addition to the biochemistry of metabolism, phototrophic growth is highly dependent on regulatory networks to coordinate growth and to relay environmental information. To this end, of particular relevance is the cyanobacterial circadian clock discovered in the late 1980s and unique among prokaryotes (Pattanayak and Rust, 2014). The cyanobacterial clock consists of an interrelated network of multifunctional components functioning in timekeeping, input and/or output mechanisms. *Synechococcus elongatus* PCC 7942 is the cyanobacterium whose clock is currently best studied. The core oscillator comprises only three proteins: KaiA, KaiB, and KaiC (Ishiura et al., 1998). KaiC exhibits an intrinsic kinase, dephosphorylation, and ATPase activity (Nishiwaki et al., 2004; Terauchi et al., 2007; Egli et al., 2012; Nishiwaki and Kondo, 2012). In complex with KaiA and KaiB, KaiC undergoes circadian Thr/Ser phosphorylation (Nakajima et al., 2005; Nishiwaki et al., 2007; Rust et al., 2007) (**Figure 6**). KaiA promotes and KaiB represses phosphorylation of KaiC (Iwasaki et al., 2002; Kitayama et al., 2003; Xu et al., 2003). The ATPase activity of KaiC is extremely weak (only 15 ATP molecules are consumed per day) and slow, determining the about 24-h period of the clock (Terauchi et al., 2007; Abe et al., 2015; Chang et al., 2015). The circadian rhythm of KaiC phosphorylation runs without transcription–translation and can even operate in a test tube (Nakajima et al., 2005; Tomita et al., 2005). *In vivo*, the KaiABC protein system works as a post-translational oscillator (PTO).

The KaiC phosphorylation rhythm has widely been studied in systems biology, and hence, a variety of mathematical models have been put forward. The first models were rather minimal to explain how sustained oscillations in phosphorylation of KaiC occur, including no intermediate steps of phosphorylation, introducing feedbacks on KaiC phosphorylation or assuming hypothetical states of KaiA, KaiB, and KaiC (Emberly and Wingreen, 2006; Kurosawa et al., 2006; Mehra et al., 2006; Axmann et al., 2007; Mori et al., 2007). Emberly and Wingreen (2006) were the first who showed theoretically that monomer shuffling between KaiC hexamers at specific clock times could explain the robustness and resilience of the circadian clock—a hypothesis stated prior to experimental evidence. Different variations of the concept, monomer shuffling, have afterward been modeled by other groups (Kageyama et al., 2006; Mehra et al., 2006; Mori et al., 2007; Yoda et al., 2007; Eguchi et al., 2008; Nagai et al., 2010). Another hypothesis of how synchrony within the Kai oscillator could be achieved stresses KaiA sequestration into KaiABC complexes (Kurosawa et al., 2006; Clodong et al., 2007; Rust et al., 2007; van Zon et al., 2007; Brettschneider et al., 2010). The consensus view has emerged that both mechanisms work in concert (Qin et al., 2010a). See **Table 5** (*in vitro*) and **Table 6** (*in vivo*) for selected models of the cyanobacterial circadian clock.

Recent studies have increased our understanding of how the core oscillator is integrated with input pathways and output pathways, which enable the clock to synchronize ("entrain") to the 24-h period of the environment and to transmit temporal information to downstream processes resulting in circadian rhythms in cellular physiology. Clock input cues involve the cellular ATP/ADP ratio, which has a direct effect on the core clock by modulating the KaiCs kinase activity (Rust et al., 2011). In particular, an increase in the ADP levels, occurring when cells are placed into darkness, inhibits further KaiC phosphorylation and thus resets the phase of oscillation to synchronize to the metabolic state of the cell. By simulating various ATP/ADP ratios that mimic different night phases, Rust et al. (2011) could recreate phase shifts in the core oscillator as seen *in vitro* and *in vivo*. For their theoretical analysis, the authors modified a previous mathematical model of the circadian clock, which was based on the rates of phosphorylation and dephosphorylation at Thr432 and Ser431 (Rust et al., 2007). In the refined model, the

physical interaction of SasA with KaiC promotes phosphotransfer to RpaA so that RpaA becomes active. During the night, the physical interaction of CikA with the KaiBC complex inhibits phosphorylation of RpaA so that RpaA becomes inactive (Gutu and OShea, 2013). Red and blue dots are phosphates at KaiC phosphorylation sites Thr432 and Ser431. Red arrows represent interactions with ATP, ADP, and oxidized quinones related to metabolic processes of phototrophic growth. The role of other input components (e.g., LdpA) and output components (e.g., LabA, RpaB) as well as the location of the quinones in the thylakoid membranes are not shown. See Table 5 (*in vitro*) and Table 6 (*in vivo*) for selected models of the cyanobacterial circadian clock.

KaiC phosphorylation. The core circadian clock generates rhythms in gene expression and cell division via the global transcriptional factor RpaA. During the day, the

KaiA-dependent kinase rates were now additionally modulated by the ratio of ATP to ATP + ADP. This model was again adapted to account for an additional ATPase activity experimentally found in the CI subunit of KaiC and required for binding of KaiB to Ser-phosphorylated KaiC (Phong et al., 2013). The KaiBC complex formation was shown to depend on an ATPase, but whose activity was insensitive to changes in the cellular ATP/ADP ratio—in contrast to the ATPase in the CII subunit of KaiC (responsible for the Thr/Ser phosphorylation reactions). The results of the combined modeling and experimental study (Phong et al., 2013) suggest that these two differently sensitive catalytic domains are responsible for the capability of the clock to receive input signals while preserving circadian rhythmicity. Depending on the specific question under investigation, both entrainment models could relatively straightforwardly be integrated with other modules of cyanobacterial physiology, with ATP as key to coupling (e.g., light reactions). Another model version already exists, which additionally accounts for the transcription and translation of clock genes as well as the feedback to the core oscillator but needs to be extended to include interactions with other cellular parameters such as the ATP/ADP ratio (Teng et al., 2013). The mathematical clock model proposed by Brettschneider et al. (2010) could equally be envisaged for the integration into larger models of phototrophic growth. The core oscillator is modeled by a larger set of ODEs (12 ODEs; for comparison, 3 ODEs in Rust et al. (2007) but includes hexamerization of KaiC, monomer shuffling, and assembly and disassembly of KaiAC, KaiBC, and KaiABC complexes. Here also, an extended model (15 ODEs) coupled to a transcription–translation circuit has been proposed (Hertel et al., 2013). Inhibition of KaiCs' kinase activity by ADP can be incorporated into the model.

#### TABLE 5 | Selected *in vitro* models of the cyanobacterial circadian clock.

*The models focus on the functioning of the post-translational oscillator (PTO). ATP, NADPH, and redox state are typically not considered as explicit variables. All models are implemented as ODEs, except otherwise noted.*

Sensing the decrease in the ATP/ADP ratio during night is assumed to allow the clock to infer the length of night (Rust et al., 2011). In two biochemical studies, Pattanayak et al. uncovered a further connection between the clock and cellular metabolism: metabolic rhythms produced by the clock (such as rhythms in glycogen abundance, which go along with changing levels of ATP/ADP) feed back to the core oscillator. These rhythms are very likely the main driving force of the clock, allowing the cells to anticipate the onset of darkness in advance (Pattanayak et al., 2014, 2015). In addition, the clock seems to be able to anticipate nightfall through the plastoquinone pool, which is part of the ETC. In particular, the plastoquinone pool embedded within the thylakoid membrane becomes transiently oxidized at the transition from day into night and binding of PQ to KaiA causing aggregation and decay that, in turn, reduces the positive effect of KaiA on KaiC phosphorylation (Wood et al., 2010; Kim et al., 2012). Other redox-sensitive input components such as CikA (circadian input kinase) and LdpA (light-dependent period

#### TABLE 6 | Selected *in vivo* models of the cyanobacterial circadian clock.

*Later models combine the transcription–translation feedback loop (TTFL) and the core PTO oscillator. As for in vitro models, ATP, NADPH, and redox state are typically not considered as explicit variables. All models are implemented as ODEs, except otherwise noted.*

A) have been identified, which reset or modulate the phase of KaiC phosphorylation cycle [reviewed by Mackey et al. (2011)]. **Figure 6** provides a schematic of possible sites of interactions. Mathematical models describing the interconnections at the molecular level have yet to be developed.

The cyanobacterial circadian clock results in genome-wide gene expression rhythms and regulates cell cycle progression relaying information via a two-component system that comprised SasA (*Synechococcus* adaptive sensor A) and RpaA (regulator of phycobilisome association A) (Liu et al., 1995; Mori et al., 1996; Takai et al., 2006; Dong et al., 2010). Rhythms of chromosome compaction and DNA topology (highly correlated with gene expression rhythms) do not hinge on SasA, pointing to the existence of other output pathways (Smith and Williams, 2006; Woelfle et al., 2007; Vijayan et al., 2009). In the positive transcriptional pathway, SasA interacts physically with KaiC and acts as a kinase toward RpaA (Takai et al., 2006; Gutu and OShea, 2013). In the course of a circadian cycle, KaiB displaces SasA from KaiC and KaiA becomes sequestered, switching KaiC into autodephosphorylation mode (**Figure 6**). The negative transcriptional output involves LabA (low amplitude and bright), CikA, and the transcriptional factor RpaB (regulator of phycobilisome associated B)—all three repressing the activity of RpaA (Taniguchi et al., 2007; Gutu and OShea, 2013; Espinosa et al., 2015). CikA, with its dual role in input and output pathways, plays a special role. As an output component, CikA competes with KaiA for binding to KaiB (Gutu and OShea, 2013; Chang et al., 2015). The binding of CikA to the KaiBC complex activates the phosphatase activity of CikA toward RpaA (Gutu and OShea, 2013) (**Figure 6**). RpaA, as both a circadian transcriptional activator and a repressor, drives global gene expression rhythms. The transcriptional output includes the regulation of clock genes, forming a transcription–translation feedback loop (TTFL) to the core oscillator (Ishiura et al., 1998). Since we are just beginning to understand the mechanistic details of the TTFL, the existing mathematical models are still highly simplified, using phenomenological assumptions as to how RpaA (Zwicker et al., 2010) or specific phospho-forms of KaiC control transcription of the *kaiBC* gene cluster (Miyoshi et al., 2007; Qin et al., 2010b; Hertel et al., 2013; Teng et al., 2013). These models reproduce the most important experimental results, although Miyoshi et al. (2007) assumed hypothetical states for KaiA, KaiB, and KaiC inconsistent with experiments. Due the relatively small numbers of variables and parameters, the ODE models of Teng et al. (2013) and Hertel et al. (2013) might be most suitable for use in larger models of phototrophic growth but require additional modifications that account for the most recent experimental findings.

## 8. REGULATION OF GENE EXPRESSION IN CYANOBACTERIA

In *Synechococcus elongatus* PCC 7942, the KaiC phosphorylation cycle targets the general transcription apparatus and thereby regulates 30–100% of the transcriptome in circadian fashion, depending on the experimental setup (Liu et al., 1995; Nakahira et al., 2004; Ito et al., 2009; Vijayan et al., 2009; Lehmann et al., 2013). The transcriptional output is regulated by multiple factors such as circadian changes in chromosomal compaction/ decompaction (Smith and Williams, 2006) involving oscillations in DNA supercoiling (Woelfle et al., 2007) as well as biochemical cascade pathways that converge to globally acting transcription factors, RpaA and RpaB and so far unknown factors (Taniguchi et al., 2007; Gutu and OShea, 2013; Paddock et al., 2013; Espinosa et al., 2015). Furthermore, it is now clear that small non-proteincoding RNAs (<200 nucleotides) play as both positive and negative regulators crucial roles in gene expression of cyanobacteria (Georg and Hess, 2011). By base-pairing with the target mRNA, small RNA molecules interfere with the ribosome binding site or other sequence stretches, and consequently alter mRNA translation and stability. This mode of regulation might explain why the proportion of cyclic proteins in diverse cyanobacteria is rather uncorrelated to that found in microarray studies (Stöckel et al., 2011; Waldbauer et al., 2012; Guerreiro et al., 2014), because not only transcriptional but also post-transcriptional (small RNAmediated) mechanisms might be active and modulate or finetune the dynamics of regulatory networks. *Synechocystis* sp. PCC 6803 possesses a large number of small non-coding RNAs, and antisense RNAs influence at least 26% of all gene transcripts in this cyanobacterium. There are several hints that the non-coding RNAs fulfill important functions in light–dark acclimation (Georg et al., 2009; Mitschke et al., 2011). A specific functional role could be clarified for some of the antisense RNAs, e.g., IsrR (Dühring et al., 2006; Legewie et al., 2008), as-flv4 (Eisenhut et al., 2012), or PsbAR2 and PsbAR3 mRNA (Sakurai et al., 2012). Yet, many identified RNA regulators still await elucidation of their functional relevance.

# 9. CONCLUSION: PUTTING THE PARTS TOGETHER

Almost all cellular functions have evolved in the context of constraints and trade-offs that can only be understood if the respective cellular and environmental context is taken into account. To this end, the construction of computational models of cellular processes not only allows us to study the inner workings of selected processes but also allows us to investigate the emergent properties that arise from interactions between these processes. The tradeoffs and interrelations within phototrophic growth are manifold: the energy required for cellular growth is derived from the photosynthetic light reactions, which themselves are a major source of reactive oxygen species (ROS) and therefore require careful balance between different electron transport pathways and alternative electron "valves". CCMs use energy and have implications for the efficiency of carbon assimilation. The cellular ATP/ADP ratio and the oxidized PQ pool relays information to the circadian clock, which affects transcriptional output and hence metabolic activity. Metabolism itself depends on cellular energy and redox potential—and must be appropriately coordinated to synthesize the right metabolites at the right time. The availability of ribosomes and amino acids, which are itself products of metabolism, affects the rates of translation of new proteins—which must be coordinated to account for damage, stability, and turnover times of proteins. In particular, the components of PSII complexes themselves are dependent on environmental conditions due to photodamage caused by ROS (Yao et al., 2012).

Many if not most of these interactions and trade-offs are still insufficiently understood. An important example is the action and the evolutionary benefit of the circadian clock. While many if not most prokaryotes live in environments with periodic diurnal cycles of light, temperature, and humidity, only cyanobacteria are known to possess a *bona fide* circadian clock. While the competitive advantage of a circadian clock in a periodic environment has been demonstrated experimentally for cyanobacteria (Woelfle et al., 2004), the precise adaptive value and the selective pressure resulting the evolution of a clock remains only partially understood. Reasoning about the possible fitness implications of a circadian clock necessarily involves considering the organisms as a whole, as exemplified in the "escape from light" hypothesis that circadian rhythmicity arose from the need to protect the organism's DNA from ultraviolet (UV) radiation, at the time unfiltered by the Earth's early atmosphere (Hut and Beersma, 2011; Lück and Westermark, 2016). A quantitative evaluation of such a hypothesis requires to contrast the energetic cost of the circadian clock with its benefits for survival and growth—a task where advanced computational models will allow for an increasingly quantitative evaluation.

While, as outlined in this contribution, a large number of computational models related to phototrophic growth are already available, also many important cellular processes are still insufficiently described. An important example is the coordination of cellular growth in response to transient darkness, starvation, or stress conditions. Only recently, iconic pathways, such as the stringent response, have been shown to be active also in cyanobacteria and to mediate a coordinated transcriptional and translational reaction to (transient) periods of darkness (Hood et al., 2016). Likewise, knowledge about the molecular and physiological mechanisms involved in the transition of cyanobacterial cells from a resting state to an active vegetative state is still incomplete (Klotz et al., 2016)—albeit crucial to understand what mechanisms causes a cyanobacterial cell to divert resources away from growth and division and toward survival until environmental conditions improve. Also, only little is known concerning the coordination of cellular metabolism with cell cycle events (Asato, 2005, 2006). In particular, cell size control and size homeostasis in bacteria is still not fully understood, with conceptual models dating back to the work of Donachie (1968). Recent work on *E. coli* and *B. subtilis* favored the "ädder" model, in which the size added between birth and division is constant for a given growth condition—as opposed to the "sizer," in which the cell actively monitors cell size, and "timer" model, in which the cell grows for a specific time before division (Taheri-Araghi et al., 2015). Corresponding studies for cyanobacteria are not yet available. Recent data have indicated that there is coupling between circadian oscillator and the cell cycle, specifically that cell cycle progression in some cyanobacteria slows during specific circadian intervals (Dong et al., 2010; Yang et al., 2010)—posing timely questions for further computational research and necessitating integrative modules of cyanobacterial growth.

Overall, there is increasing interest in whole-cell models to understand cellular trade-offs and functions in the context of a living cells. The construction of integrative computational models to predict phenotype from genotype has gained momentum with a first whole-cell model of the life cycle of the human pathogen *Mycoplasma genitalium*—based on a subdivision of cell functionality into modules (Karr et al., 2012). A similar effort was already undertaken for cyanobacteria to explain fitness advantage conveyed by a circadian clock (Hellweger, 2010)—an approach that can be regarded as a prototype for the path outlined in this contribution.

Distinct from other approaches to whole-cell models, however, we argue that it is unlikely that a single universal model—a model that spans all scales from intracellular to intercellular to properties of ecosystems—will fulfill all requirements needed to describe cellular growth. Rather, we envision a modular approach. Depending on the research question, different temporal and spatial scales must be considered. The task is then to derive an appropriate representation of cellular processes that accounts for the spatial and temporal scales involved. The derived submodels should be consistent with more fine-grained representation and their construction must be informed by knowledge how processes interact and which trade-offs are relevant for a specific research question. We note that, besides the biological challenges, such a strategy also entails major computational challenges. As yet, the annotation of computational models is often poor. That is, the biochemical identity of model variables is not defined in a computer-readable format, and hence, merging of models typically requires extensive manual curation (Krause et al., 2010). While standardized exchange formats for computational models, such as the Systems Biology Markup Language, SBML (Hucka et al., 2003), are available for more than a decade, they are not commonly applied outside the Systems Biology community. For example, as yet, none of the models of the cyanobacterial circadian clock are available from the BioModels database—a major resource for computational models of biological processes (Li et al., 2010). In addition, only few computational tools allow for the integration of different modeling concepts, such as constraintbased and kinetic models.

Notwithstanding these challenges in computational methodology, we expect a growing library of models related to cyanobacterial growth that can inform and guide research how phototrophic organisms, such as cyanobacteria, adapt to complex environments. These models must increasingly adhere to common standards and should be made available on open platforms, such as the BioModels database (Li et al., 2010) and e-cyanobacterium.org (Klement et al., 2013). Given the recent

# REFERENCES

advances in model development and annotation, computational modeling will undoubtedly play a key role in understanding trade-offs and adaptations in cyanobacteria. In the beginning, integrated models of cyanobacterial growth will be still idealistic, crude, and most certainly incomplete. But, again referring to Neidhardt (1999), "it is only through such modeling of wholesystem behavior—that is, of growth—that one will learn how near and how far our knowledge takes us toward understanding the living cell."

# AUTHOR CONTRIBUTIONS

SW and RS conducted the research and wrote the article. Both authors approved the final manuscript.

# ACKNOWLEDGMENTS

We would like to thank Jiří Jablonský for discussing his work on modeling phtosynthetic light reactions and carbon metabolism, as well as Robert L. Burnap for comments and sharing insight about cyanobacterial physiology and growth.

# FUNDING

This work is a contribution from the research grant *CyanoGrowth* funded by the German Federal Ministry of Education and Research as part of the e:Bio Innovationswettbewerb Systembiologie [e:Bio systems biology innovation competition] initiative (reference: FKZ 0316192).

should we take into account photo-acclimation and variable stoichiometry in oligotrophic areas? *J. Mar. Syst.* 125, 29–40. doi:10.1016/j.jmarsys.2012.12.010

and in the isolated reaction center: pheophytin is the primary electron acceptor. *Proc. Natl. Acad. Sci. U.S.A.* 103, 6895–6900. doi:10.1073/pnas.0505371103

photoautotrophic and heterotrophic growth conditions. *Photosynth. Res.* 118, 167–179. doi:10.1007/s11120-013-9943-x

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Copyright © 2016 Westermark and Steuer. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

**146**

# APPENDIX

The differential equations that define the model shown in **Figure 4** are d[PSU\*]/d*t* = *v*<sup>1</sup> − *v*<sup>2</sup> − *vd* and d[PSU*d*]/d*t* = *vd* − *vr*, together with the conservation relationship [PSU] + [PSU\*] + [PSU*d*] = [PSUtotal]. The reactions rates describe the activation of the photosynthetic unit (PSU), *v*<sup>1</sup> = *k*1⋅*I*⋅[PSU], the return to the open state, *v*<sup>2</sup> = *k*2⋅[PSU\*], the transition to a photodamaged state, *vd* = *kd*⋅*I*⋅[PSU\*], and the recovery, *vr* = *kr*⋅[PSU*d*]. The rate of activation, *v*1, and photodamage, *vd*, both depend on the light intensity, *I* (with *kd* ≪ *k*1). Parameters used are *k*<sup>1</sup> = 10.0, *k*<sup>2</sup> = 10.0, *kr* = 1.0, [PSUtotal] = 1.0, and *kd* as indicated in **Figure 4**. The rate of oxygen evolution (*v*2) obtained from an (analytical) steady-state solution of the ODEs is shown. Parameters are arbitrary and do not reflect actual values found in cyanobacteria.

# Customized Steady-State Constraints for Parameter Estimation in Non-Linear Ordinary Differential Equation Models

#### Marcus Rosenblatt <sup>1</sup> \*, Jens Timmer 1, 2, 3 and Daniel Kaschek <sup>1</sup>

1 Institute of Physics, Albert Ludwig University of Freiburg, Freiburg, Germany, <sup>2</sup> Freiburg Centre for Systems Biology, Albert Ludwig University of Freiburg, Freiburg, Germany, <sup>3</sup> BIOSS Centre for Biological Signaling Studies, Albert Ludwig University of Freiburg, Freiburg, Germany

Ordinary differential equation models have become a wide-spread approach to analyze dynamical systems and understand underlying mechanisms. Model parameters are often unknown and have to be estimated from experimental data, e.g., by maximum-likelihood estimation. In particular, models of biological systems contain a large number of parameters. To reduce the dimensionality of the parameter space, steady-state information is incorporated in the parameter estimation process. For non-linear models, analytical steady-state calculation typically leads to higher-order polynomial equations for which no closed-form solutions can be obtained. This can be circumvented by solving the steady-state equations for kinetic parameters, which results in a linear equation system with comparatively simple solutions. At the same time multiplicity of steady-state solutions is avoided, which otherwise is problematic for optimization. When solved for kinetic parameters, however, steady-state constraints tend to become negative for particular model specifications, thus, generating new types of optimization problems. Here, we present an algorithm based on graph theory that derives non-negative, analytical steady-state expressions by stepwise removal of cyclic dependencies between dynamical variables. The algorithm avoids multiple steady-state solutions by construction. We show that our method is applicable to most common classes of biochemical reaction networks containing inhibition terms, mass-action and Hill-type kinetic equations. Comparing the performance of parameter estimation for different analytical and numerical methods of incorporating steady-state information, we show that our approach is especially well-tailored to guarantee a high success rate of optimization.

Keywords: non-linear ODE models, parameter estimation, biochemical reaction networks, steady-state, positive solutions, multiplicity, multi-stability, success rate

# 1. INTRODUCTION

Dynamical systems are frequently modeled by systems of ordinary differential equations (ODEs). Homogeneously distributed molecules are treated as continuous quantities interacting with each other according to kinetic laws, e.g., mass-action or Michaelis-Menten kinetics.

#### Edited by:

Julio Vera González, University Hospital Erlangen, Germany

#### Reviewed by:

Irene Otero Muras, Consejo Superior de Investigaciones Científicas, Spain Luis L. Fonseca, Georgia Institute of Technology, USA

\*Correspondence:

Marcus Rosenblatt<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Systems Biology, a section of the journal Frontiers in Cell and Developmental Biology

> Received: 04 December 2015 Accepted: 21 April 2016 Published: 11 May 2016

#### Citation:

Rosenblatt M, Timmer J and Kaschek D (2016) Customized Steady-State Constraints for Parameter Estimation in Non-Linear Ordinary Differential Equation Models. Front. Cell Dev. Biol. 4:41. doi: 10.3389/fcell.2016.00041 A typical ODE system

$$
\dot{\boldsymbol{x}} = \boldsymbol{f}(\boldsymbol{x}, \boldsymbol{p}, \boldsymbol{u}(t)) \,, \qquad \boldsymbol{\kappa}(0) = \boldsymbol{x}\_0 \tag{1}
$$

determines the time-evolution of an N-dimensional state vector x(t). Here, p ∈ R M <sup>+</sup> denotes the M-dimensional vector of nonnegative kinetic parameters. The vector x<sup>0</sup> ∈ R N +,0 , where R+,<sup>0</sup> = R<sup>+</sup> ∪ {0}, gives the set of initial conditions. The kinetic parameters and initial conditions together span the space of model parameters θ = (p, x0). The explicit time-dependency via u(t) corresponds to external driving forces, like drug stimuli in biological dynamic systems.

In many fields where ODE models are used, parameter values are not a priori known and have to be estimated from experimental data. Commonly, this is achieved by minimizing an objective function g(θ , D) that penalizes weighted differences between model prediction x(t) and data D, e.g., by maximumlikelihood estimation. For the case of non-linear ODE systems, several local optima may exist. In order to find the global optimum, several optimization methods, e.g., particle swarm optimizers (Peng et al., 2010) or simulated annealing (Xiang and Gong, 2000), include stochasticity to escape local minima. Compared to that, deterministic algorithms may stick to local optima during optimization. On the other hand, gradient and Hessian information of the objective function can be incorporated, increasing the performance of optimization by a multiple. Combining the advantages of derivative-based optimization and random sampling, a multi-start deterministic optimization approach has proven to yield superior overall performance for our problem class (Raue et al., 2013). Throughout this work, we perform optimization by means of a trust-region optimizer from multiple starting positions.

Specially in models of biological systems, available data is sparse and parameters are often non-identifiable. Apart from that, the high-dimensional parameter space hampers parameter sampling. In order to reduce the number of parameters, the system is assumed to initially (t = 0) be in a steady-state which is determined by the constraint equation

$$f(\mathbf{x}\_0, p, \mathbf{0}) = \mathbf{0}.\tag{2}$$

As a standard approach, the steady-state constraint is solved for the initial values x0. Since Equation (2) is in general non-linear, this may lead to higher-order polynomial equations for which no general solution is available. Even for a rather simple case of quadratic or cubic equations, solutions are not unique and optimization would have to be performed for all possibilities. Another aspect of steady-state calculation are negative solutions for x<sup>0</sup> and p that appear for certain model specifications. Negative solutions are not only contradicting the biological setting with positively defined concentrations and kinetic parameters but also constitute a problem for optimization. Negative parameter values change the sign of damping terms of the ODE's right-hand side which might lead to rapidly growing solutions and an abort of the optimization before an optimum was reached.

In order to obtain a high convergence probability for the optimization of randomly chosen initial parameter samples, our aim is to derive non-negative, analytical steady-state expressions, while multiple steady-state expressions are likewise circumvented by a proper choice of kinetic and initial value parameters for which Equation (2) is solved.

Over the last decades, steady-state analysis has been addressed by many algorithms and methods. In the following, we give an overview of existing approaches and summarize their applicability to different types of model equations with a special focus on parameter estimation in ODE models, see **Table 1**.

The earliest-proposed algorithm for deriving steady-states in enzyme-catalyzed systems being described by simple mass-action rules was developed by King and Altman (1956). In the original paper, however, interactions that do not involve the enzyme were not allowed which prohibits applicability to most of today's systems with proteins mediating the activation of other proteins without being part of the reaction. After Chemical Reaction Network Theory (CRNT) was formulated (Horn and Jackson, 1972; Feinberg, 1979), the method of King and Altman has been improved by graph theory (Chou, 1990) and extended to special subclasses of CRNs, e.g., layered signaling cascades (Feliu et al., 2012) and post-translational modification networks (Feliu and Wiuf, 2013). The same authors also published a more general approach for CRNs in Feliu and Wiuf (2012). Here, a set of core variables is introduced serving for a parametrization of the steady-states whereby non-negative solutions are guaranteed due to graph-theoretical arguments.

Another approach developed by Halasz et al. (2013), introduces bilinearities of the system as new variables leading to a linearized system solvable by application of Cramer's rule. The number of bilinearities, however, is restricted and negative steady-state solutions are not prevented.

All mentioned approaches deal with steady-state analysis for CRNs based on mass-action rules. However, modern modeling approaches often make use of special reaction types such as inhibition, Michaelis-Menten or Hill kinetics that cannot be included into standard CRNT without changing the model structure and introducing new dynamical variables. In the approach of Halasz et al. (2013), inhibition and Michaelis-Menten terms can easily be integrated by multiplying the corresponding steady-state equation by the denominator of the rate expression. However, since a state variable is contained in the denominator, this can increase the number of bilinearities significantly.

In order to avoid problems of higher-order polynomial equations, steady-state equations can be solved not only for initial value- but also for kinetic parameters, which is done in the steady-state solver py-substitution developed by Loriaux et al. (2013). From N initial values and M kinetic parameters, a set of N variables is chosen that have to be fixed in Equation (2). In doing so, a lot of freedom is incorporated into the solution. In fact, py-substitution is able to solve the very most steady-state equation systems, since in principle N kinetic parameters could be chosen as fixed variables directly leading to a simple linear equation system.

Complementary to analytical approaches, steady-state information can be incorporated into the system by numerically computing the initial conditions during each optimization step. Even gradient information that is necessary for efficient optimization is available by means of the implicit function theorem. A numerical incorporation of steady-state information has the advantage that the complexity of the underlying equation system is in principal not restricted. Furthermore, the implementation remains untouched when model equations are changed. However, convergence of the numerical steady-state calculation is not guaranteed and issues of multiple steady-states cannot be controlled.

In the following section, we present a method to derive nonnegative steady-state expressions for a large class of nonlinear ODE models that are based on biochemical reactions. Our approach picks up the idea of solving for kinetic parameters in order to derive unique and simple steady-state expressions. Due to the structure of the ODE system, solving for kinetic parameters often leads to potentially negative steady-state solutions, depending on the point of evaluation in parameter space. By introducing appropriate parameter transformations and exploiting the given model structure, our approach guarantees a non-negative solution space. In the Results section, we show how different steady-state parameterizations influence the optimization procedure and compare our approach to the standard approach of solving for initial value parameters as well as to a numerical steady-state approach.

# 2. METHODS

### 2.1. Theoretical Background

Let us consider a model f as an N-dimensional ODE system x˙ = f(x, p) with states x, parameters p and no external driving forces, i.e., the ODE is autonomous. We write f as a matrix product

$$f(\mathbf{x}, p) = \mathbf{S} \cdot F(\mathbf{x}, p) \,, \tag{3}$$

of the N × M-dimensional stoichiometry matrix S and the M-dimensional flux vector F which depends on states and parameters. For the entries of the flux vector, we allow rational functions of x and p including e.g., mass-action, inhibition, Michaelis-Menten and Hill-Type kinetics. **Table 2** gives an overview of the main reaction types covered by the presented steady-state approach.

We assume that each single flux F<sup>l</sup> is proportional to some flux parameter k<sup>l</sup> and can be written as

$$F\_l = k\_l \cdot G\_l(\mathbf{x}, q) \,, \tag{4}$$

where the function G<sup>l</sup> only depends on the states and a set of additional parameters q taken from the set of all model parameters p. The union of flux parameters k and additional model parameters q coincides with the parameter set p. Typically all reaction types described by CRNT only need one flux parameter and do not contribute to q, however, inhibition terms and Michaelis-Menten kinetics contain at least one additional parameter and Hill kinetics even two.

The signs of the entries of the stoichiometry matrix S determine whether a flux contributes as an in- or an outflux to the time evolution of the corresponding state. We assume that each outflux is at least linearly dependent on the corresponding

All types are covered by our steady-state approach.

state, as being always the case for mass-action systems. By means of Equation (2), each initial value x0,<sup>i</sup> ≡ x<sup>i</sup> is directly related with a steady-state equation of the form

$$0 = \sum \mathbf{in}\_i - \mathbf{x}\_i \cdot \sum \mathbf{out}\_i \,, \tag{5}$$

where in<sup>i</sup> and out<sup>i</sup> constitute functions of states and parameters. For a majority of reaction types used in ODE models, the fluxes in<sup>i</sup> and out<sup>i</sup> are independent of x<sup>i</sup> , compare **Table 2**. In these cases, Equation (5) is linear in x<sup>i</sup> and has the solution x<sup>i</sup> = P P ini out<sup>i</sup> . However, if the fluxes in<sup>i</sup> or out<sup>i</sup> still depend on x<sup>i</sup> , e.g., reaction 12 in **Table 2** for the case of self-activation or reaction 3 with an outflux being quadratic in x<sup>i</sup> , Equation (5) might be non-linear in x<sup>i</sup> .

In order to solve the complete set of steady-state equations, we analyze their specific structure by means of graph theory. We therefore rewrite Equation (5) as

$$\mathbf{x}\_{i} = \frac{\sum \mathbf{in}\_{i}}{\sum \mathbf{out}\_{i}} \tag{6}$$

and summarize appearances of states on the right-hand side of Equation (6). Here, the set of states is defined by the set of dynamic variables x that we want to fix by the steady-state determination. Once a dynamic variable is fixed by a nonnegative expression or treated as a free parameter, it is removed from the set of states.

**Definition 1:** A head of state x<sup>i</sup> is a state x<sup>j</sup> that appears on the right-hand side of Equation (6). By h(xi), we refer to the set of heads for a specific state x<sup>i</sup> . In particular, x<sup>i</sup> can itself be a part of h(xi).

**Proposition 1:** If non-negative steady-state solutions for all heads of x<sup>i</sup> are known, a non-negative steady-state solution for x<sup>i</sup> can directly be obtained by Equation (6). This holds especially, if the set h(xi) is empty.

**Definition 2:** The adjacency matrix M(f) of an ODE model f(x, p) with states x and parameters p is an N × N matrix with entries

$$M\_{ji} = \begin{cases} 1, \text{if} & \mathbf{x}\_j \in h(\mathbf{x}\_i), \\ \mathbf{0}, \text{ else}. \end{cases}$$

Each dM-dimensional adjacency matrix M defines a directed graph G<sup>M</sup> with nodes x<sup>1</sup> to xd<sup>M</sup> which we call steady-state graph. Each non-zero entry of M corresponds to a directed edge(xj, xi) implying that x<sup>j</sup> occurs in the steady-state expression of x<sup>i</sup> , i.e., Equation (6). A non-zero diagonal entry Mii reflects that the corresponding steady-state equation is non-linear in x<sup>i</sup> .

#### 2.2. Splitting Cycles

The specific structure of the steady-state graph enables to solve the steady-state equations step-by-step as is shown in the following.

**Definition 3:** A cycle of a steady-state graph is a path through the graph along its edges with equal starting and end point. Here, we allow cycles of length one arising from non-zero diagonal entries in the adjacency matrix M.

**Definition 4:** Graphs that do not contain cycles are called tree-like.

**Proposition 2:** If a steady-state graph of an N-dimensional model f is tree-like, non-negative steady-state solutions can be obtained for all x<sup>i</sup> inside the graph.

Proof: For any tree-like steady-state graph, there exists at least one root, i.e., state without head, called x<sup>r</sup> . Since h(xr) = ∅, the corresponding steady-state expression can be obtained by Proposition 1. In doing so, x<sup>r</sup> is removed from the steady-state graph and a new state serves as root for which Proposition 1 again gives the corresponding steady-state expression. By iteratively applying Proposition 1 for each of the N nodes, the complete steady-state solution is obtained.

Considering Proposition 2, it is clear that solving the steadystate constraint Equation (2) for the set of initial values only becomes intricate, if there are cycles inside the steady-state graph such that higher-order polynomial equations arise. The idea of our steady-state approach is to split all these cycles step-by-step such that Proposition 2 can ultimately be applied to the remaining graph.

The simplest way of splitting a cycle is by means of a conserved quantity (CQ) of the system arising from the stoichiometry. A general introduction can be found in Loriaux et al. (2013) or Halasz et al. (2013). The following definition restricts to the properties being relevant for the presented approach.

**Definition 5:** A conserved quantity (CQ) of the model f is an expression of states and parameters which remains constant during the time-evolution of f . For each CQ, the number of independent steady-state equations is reduced by one implying that one state or flux parameter that appears in the CQ can be chosen freely. If all CQs can be derived from the stoichiometry matrix, the number of CQs is given by ncq = N − RS, with the model size N and the rank of the stoichiometry matrix RS. The cases for which ncq > N − R<sup>S</sup> are discussed in Section 2.5.

In order to split a cycle by a CQ, one of the states, x<sup>c</sup> , that appears both inside the cycle and in the CQ is chosen freely. The corresponding steady-state equation is removed from Equation (2), whereby the number of independent steady-state equations remains constant. Since the state x<sup>c</sup> is treated as a free variable, all edges originating from and leading to x<sup>c</sup> can be removed from the steady-state graph and the considered cycle is split. Note, that each CQ can only be used once.

If no states inside the cycle appear in CQs, the cycle can be split by solving the steady-state equation of a specific cycle state xi for a flux parameter k<sup>l</sup> . By means of Equations (3) and (4), the steady-state expression of k<sup>l</sup> holds

$$k\_l = \frac{-1}{S\_{il}G\_l(\mathbf{x}, q)} \sum\_{j \neq l} S\_{lj} k\_j G\_l(\mathbf{x}, q) \,. \tag{7}$$

**Proposition 3:** Let nk<sup>l</sup> be the number of appearances of the flux parameter k<sup>l</sup> , see Equation (4), inside the steady-state constraint Equation (2). Then nk<sup>l</sup> coincides with the number of non-zero entries in the l-th column of the stoichiometry matrix S.

Unless k<sup>l</sup> does not appear in other steady-state equations, i.e., nk<sup>l</sup> = 1, the considered cycle is removed from the steady-state graph without affecting other parts of the graph. However, if nk<sup>l</sup> > 1, all further appearances have to be substituted by Equation (7) which creates new edges inside the steady-state graph and possibly even new cycles. In order to keep the structure as simple as possible, flux parameters with nk<sup>l</sup> = 1 play a special role.

#### 2.3. Enforcing Positivity

Although solving for flux parameters implies linear equations and therefore structurally simple steady-state expressions, the solutions are often negative for certain model specifications. Here, we show how positivity of the expressions can be guaranteed by appropriate transformations.

The steady-state expression, Equation (7), of the flux parameter k<sup>l</sup> was derived by means of the steady-state equation of x<sup>i</sup> . The expression contains minus signs if and only if at least one of the stoichiometry entries Si,j6=<sup>l</sup> has the same sign as Sil, In this case, at least one further flux contributes to x<sup>i</sup> with the same sign as F<sup>l</sup> = klG<sup>l</sup> , namely as an in- or outflux.

**Definition 6:** For the steady-state equation of x<sup>i</sup> , Equation (5), we define µ<sup>i</sup> and ν<sup>i</sup> as the number of in- and outfluxes, respectively. Furthermore, we define the dimension of the state x<sup>i</sup> as the minimum dim(xi) = min (µi, νi).

If dim(xi) = 1, a non-negative steady-state expression is obtained by solving for the particular flux parameter being the only in- or outflux, compare first examples in **Table 3**.

If dim(xi) > 1, positivity can be enforced by performing an appropriate parameter transformation. In order to do so, we divide the fluxes contributing to x<sup>i</sup> into influxes Fin,<sup>1</sup> . . . Fin,µ<sup>i</sup> and outfluxes Fout,<sup>1</sup> . . . Fout,ν<sup>i</sup> . Then, Equation (5) reads

$$0 = \sum\_{j=1}^{\mu\_l} F\_{\text{in},j} - \sum\_{l=1}^{\upsilon\_l} F\_{\text{out},l} \cdot \tag{8}$$

Let us assume that we want to solve Equation (8) for the influx parameter kin,<sup>1</sup> = F in,<sup>1</sup> <sup>G</sup>in,<sup>1</sup> . We perform a variable transformation by defining the ratio between the remaining influxes and Fin,<sup>1</sup> as

$$r\_z = \frac{F\_{\text{in},z}}{F\_{\text{in},1}} = k\_{\text{in},z} \cdot \frac{G\_{\text{in},z}}{k\_{\text{in},1}G\_{\text{in},1}} \qquad \text{for} \quad z = 2, \ldots, \mu\_i, \tag{9}$$

TABLE 3 | Examples of solving steady-state equations for flux parameters.

where the r<sup>z</sup> replace the kinetic parameters k<sup>2</sup> to kµ<sup>i</sup> . By means of Equations (8) and (9), we obtain

$$\sum\_{l} F\_{\text{out},l} = F\_{\text{in},1} \cdot \left(1 + \sum\_{j=2}^{\mu\_i} r\_{\text{in},j}\right)^2$$

and therefore

$$k\_{\rm in,1} = \frac{1}{G\_{\rm in,1}} \sum\_{l} F\_{\rm out,l} \cdot \frac{1}{1 + \sum\_{j=2}^{\mu\_i} r\_j}. \tag{10}$$

Since Gin,<sup>1</sup> and Fout,<sup>l</sup> are positive and Equation (10) is a sum of positive contributions, a non-negative steady-state expression for kin,<sup>1</sup> is guaranteed. By means of Equation (9), the remaining flux parameters have to be substituted by the non-negative expressions

$$k\_{\rm in,z} = \frac{1}{G\_{\rm in,z}} \sum\_{l} F\_{\rm out,l} \cdot \frac{r\_{\rm z}}{1 + \sum\_{j=2}^{\mu\_i} r\_j} \qquad \text{for} \quad z = 2, \ldots, \mu\_i. \tag{11}$$

For an outflux parameter, we analogously obtain

$$k\_{\text{out},1} = \frac{1}{G\_{\text{out},1}} \sum\_{j} F\_{\text{in},j} \cdot \frac{1}{1 + \sum\_{l=1}^{\nu\_l} r\_l} \qquad \text{and} \tag{12}$$

$$k\_{\text{out},z} = \frac{1}{G\_{\text{out},z}} \sum\_{j} F\_{\text{in},j} \cdot \frac{r\_z}{1 + \sum\_{l=1}^{\upsilon\_l} r\_l} \qquad \text{for} \quad z = 2, \dots, \upsilon\_l.$$

(13)

### 2.4. Algorithm for Steady-State Determination

In the previous sections, we showed how simple steady-state expressions can be obtained (Section 2.2), while positivity is likewise guaranteed (Section 2.3). In order to split one cycle of the steady-state graph and solve for a flux parameter, a pair (xi, kj) of state and flux parameter has to be chosen, which is not unique. In the following, we suggest an algorithm based on a classification of such pairs.

According to Definitions 5, 6 and Proposition 3, we associate each pair with one of four different types:

$$(\boldsymbol{x}\_i, \boldsymbol{k}\_j) \equiv \begin{cases} \text{Type 0, } \text{ if } \boldsymbol{x}\_i \text{ appears in a CQ} \\ \text{Type 1, } \boldsymbol{n}\_{k\_j} = 1 \quad \text{and} \quad \dim(\boldsymbol{x}\_i) = 1 \\ \text{Type 2, } \boldsymbol{n}\_{k\_j} = 1 \quad \text{and} \quad \dim(\boldsymbol{x}\_i) > 1 \\ \text{Type 3, else } . \end{cases}$$

**Figure 1** shows a flowchart of the algorithm. At first, the set of CQs is computed for the ODE system serving as an input for the algorithm. If the graph is tree-like, the remaining equations are obtained according to Proposition 2 and the complete set of steady-state equations is returned.

In the case of a pair of Type 0, the cycle can simply be removed by interpreting the corresponding state as a free variable. The CQ that is thereby used is removed from the set of CQs and cannot

further contribute to the steady-state determination. Here, the flux parameters remain unaffected.

Multiplication of Equation (14) by an N × N-matrix M yields

$$\stackrel{\bullet}{\text{Unless the cycle cannot be directly split by means of a CQ, the corresponding steady-state equation, Equation (5), is solved for one of the }\mu\_i + \nu\_i\text{ flux parameters by use of a pair of Type 1,2 or 3. In order to keep the steady-state solution as simple as possible, pairs of Type 1 are preferred, since this enables to split the cycle both without substituting the flux parameter by its steady-state expression, Equation (7), and without introducing flux ratios as new parameters, Equation (9).$$

If no pairs of Type 1 are available, the algorithm scans the steady-state graph for pairs of Type 2. In this case a parameter transformation is necessary in order to guarantee positivity of the solution. However, the flux parameter does not appear any more in the system and therefore has not to be substituted. In all three cases, Type 0, 1, or 2, the number of cycles of the steady-state graph is reduced.

If pairs of Type 2 are also not available, all pairs are of Type 3. In this case, it is not a priori clear which pair is the best choice. As a simply revisable choice, the algorithm then solves the steady-state equation of the state with minimal dimension. Subsequently, all further appearances of flux parameters have to be replaced by their particular transformation, Equations (11) or (13).

# 2.5. Calculating the Conserved Quantities and Simplifying the Stoichiometry Matrix

In order to find CQs of the ODE system, linear combinations of rows of the stoichiometry matrix S can be analyzed. According to Equation (3), the N-dimensional ODE system can be written as

$$
\dot{\mathfrak{x}} = \mathbb{S} \cdot F(\mathfrak{x}, \mathfrak{p}) \,. \tag{14}
$$

$$\mathcal{M} \cdot \dot{\mathfrak{x}} = \tilde{\mathcal{S}} \cdot F(\mathfrak{x}, \mathfrak{p}) \,, \tag{15}$$

where the matrix S˜ = M · S defines linear combinations of rows of S. For each row S˜ <sup>i</sup> that is equal to zero, the quantity Mi · x = P <sup>j</sup> Mijx<sup>j</sup> is conserved.

In fact, each set of linearly dependent rows of S implies a CQ. For some ODE systems, however, not all CQs can be derived from S without accounting for the flux vector F. Equation (14) can be written

$$
\dot{\mathfrak{x}} = \mathbf{C}(\mathfrak{p}, \mathfrak{x}) \cdot \mathfrak{x} \,. \tag{16}
$$

where C(p, x) is an N × N-matrix dependent on the parameters and states. Analyzing linear dependencies of C, all CQs of the form

$$\sum\_{j} a\_{j}(\mathbf{x}, \boldsymbol{\rho}) \cdot \mathbf{x}\_{j} = \text{const.} \tag{17}$$

can be found, where the coefficients a<sup>j</sup> might depend on states and parameters.

In order to determine symbolic expressions for the a<sup>j</sup> , we transpose the matrix C and numerically search for linearly dependent columns. All parameters and states appearing in C T are replaced by random values to obtain a numeric matrix C T ran for which a QR-decomposition is performed. The matrix R constitutes an upper triangular matrix, where the number of nonempty rows corresponds to the rank of Cran. The first column R<sup>ℓ</sup> with Rℓℓ = 0 is a linear combination of the columns Rj<ℓ. Therefore, also the column Cℓ is a linear combination of the columns Cj<ℓ implying that the equation system P j ajC<sup>j</sup> = 0 has a solution for the a<sup>j</sup> with aj>ℓ = 0 which can be calculated symbolically. Thus, the quantity a T · x is conserved. Once a CQ has been found, one of the corresponding linearly dependent rows of the stoichiometry matrix is removed and the procedure is repeated. In most ODE systems, all CQs of the system can be obtained in that way. For all other cases, our Python code provides a possibility to manually specify CQs.

The idea of taking linear combinations of the stoichiometry matrix S, see Equations (14) and (15), can be augmented to simplify S for the calculation of steady-state expressions. For each matrix M, the original steady-state constraint, S · F = 0, is replaced by a new set of steady-state equations, S˜ · F = 0. With a clever choice of M, these new steady-state equations might be structurally simpler than the original ones. With respect to our proposed algorithm, the matrix M should (1) minimize the overall number of entries in the new stoichiometry matrix S˜ and (2) prevent the creation of new cycles. In practice, the idea of linearly combining rows of the stoichiometry matrix can lead to structurally simpler steady-state expressions as we show by means of a small example in the Supplementary Material.

## 2.6. Numerically Computed Steady-States

Besides calculating steady-states analytically, roots of the steadystate constraint, Equation (2), can be computed numerically during each step of the optimization. Here, we perform Newton's method which is fast compared to the time of the ODE integration. The gradient information that is necessary within our deterministic optimization scheme is determined by the implicit function theorem, i.e., given the steady-state constraint

$$f(\mathbf{x}\_0(p), p) = 0,$$

we derive the equation with respect to p and obtain

$$0 = \frac{\partial f}{\partial \mathbf{x}\_0} \cdot \frac{\partial \mathbf{x}\_0}{\partial p} + \frac{\partial f}{\partial p} \qquad \Longrightarrow \qquad \frac{\partial \mathbf{x}\_0}{\partial p} = -\left(\frac{\partial f}{\partial \mathbf{x}\_0}\right)^{-1} \frac{\partial f}{\partial p} \cdot \mathbf{x}\_0$$

### 2.7. Technical Remarks

The steady-state algorithm was implemented in Python by use of the libraries numpy and sympy. It can either be downloaded from the author's homepage as a Python code or can be used from within the R-packages dMod/cOde available from https://github.com/dkaschek/. Simulation of data and parameter estimation with analytical and numerical steady-states were performed in dmod.

# 3. RESULTS

When calculating steady-state expressions for parameter estimation of ODE systems, several aspects have to be considered simultaneously. Most importantly, the parameter space is to be reduced as far as possible. Therefore, all available steady-state constraints should be taken into account. Since solving for state variables often leads to higher-order equations for which solutions are difficult to obtain, one has at least partially to solve for kinetic parameters. In doing so, the steady-state expressions often lead to negative parameter values for certain model specifications.

Due to mass balance, outfluxes contribute with a minus sign to the time derivative of the corresponding state. Provided that outflux rates are proportional to positive powers of their states, they contribute damping terms to the time-evolution of the state. However, if for a certain model specification the corresponding flux parameter is negative, the sign of the outflux term becomes positive which leads to an exploding model trajectory for the state.

In Section 3.1, we show how our steady-state approach determines simple steady-state equations for systems that lead to higher-order equations when solved for the state variables. In Section 3.2, we show how steady-state expressions with negative realizations lead to optimization problems and a significantly lower success rate, i.e., the probability to converge to a local or the global optimum. Non-linear ODE systems often have several steady-state solutions, when the steady-state equations are solved for the state variables. For parameter estimation, multiple steadystates constitute a problem, since all possible realizations have in principle to be followed up. By solving for kinetic parameters, our steady-state approach likewise avoids multiple solutions and improves the optimization as we show in Section 3.3.

# 3.1. Determination of Non-Negative Steady-State Expressions

To show the applicability of the presented steady-state approach, we investigate a toy model with six state variables and nine reactions of the form

∅ <sup>k</sup><sup>0</sup> −→ A A <sup>k</sup><sup>1</sup> −→ B B <sup>k</sup><sup>2</sup> −→ <sup>A</sup> A + A <sup>k</sup><sup>3</sup> −→ C C <sup>k</sup><sup>4</sup> −→ ∅ <sup>B</sup> <sup>+</sup> <sup>C</sup> <sup>k</sup><sup>5</sup> −→ <sup>D</sup> D + G <sup>k</sup><sup>6</sup> −→ F B <sup>k</sup>7·<sup>F</sup> −−→ ∅ F <sup>k</sup><sup>8</sup> −→ <sup>G</sup>.

All reactions satisfy the law of mass action, the degradation of B is mediated by F. With these assumptions, one obtains the following ODE system

$$\begin{aligned} \dot{A} &= k\_0 + k\_2 B - k\_1 A - k\_3 A^2 \\ \dot{B} &= k\_1 A - k\_2 B - k\_5 B C - k\_7 B F \\ \dot{C} &= k\_3 A^2 - k\_5 B C - k\_4 C \\ \dot{D} &= k\_5 B C - k\_6 D G \\ \dot{F} &= k\_6 D G - k\_8 F \\ \dot{G} &= k\_8 F - k\_6 D G. \end{aligned}$$

The system contains one conserved quantity F + G = const, reflecting that the steady-state equations of F and G are not independent from each other. Therefore, the number of variables that have to be fixed by the steady-state is five. In order to obtain the corresponding steady-state equations, all time-derivatives of the states are set to zero. Although the single equations of this system are of degree two or lower, solving for the states leads to a sixth order polynomial equation, see Supplementary Material, for which no closed-form solution is available.

TABLE 4 | Steps of steady-state determination for a model with six states and eight reactions.

**Table 4** summarizes how our steady-state solver determines a non-negative steady-state solution by partially solving for flux parameters. During the first loop, the cycle [A, A] of state A to itself is split. The pairs of A and its contributing flux parameters are all of Type 3, since there are two influx- and two outflux parameters of which at least one is appearing in the other steadystate equations, e.g., in the equation of B. The equation of A is solved for the influx parameter k0, whereby k<sup>2</sup> is transformed and replaced by the new free parameter r<sup>1</sup> = k2B/k0, see first loop in **Table 4**. The appearance of k<sup>2</sup> in the equation of B is substituted, whereas k<sup>0</sup> has no further appearances.

Whereas the state A is removed from the steady-state graph, the state B has become a new head of itself in consequence of the substitution. In the second loop, this new cycle [B, B] is split by solving the equation of B for the flux parameter k1. Here, the pair (B, k1) is of Type 1, since k<sup>1</sup> is the only influx parameter and not appearing in the remaining equations.

In the next loop, the algorithm splits the cycle [D, G, D] by taking G as a free parameter, since it is part of the conserved quantity F + G. The remaining steady-state graph in the last loop is tree-like and therefore the steady-state equations can be derived according to Proposition 2 starting with C which in this case serves as the root of the graph.

For simplification of writing, our steady-state solver outputs the equations in a specific order where fixed states or parameters may still appear in the equations below. In order to obtain a complete independent set of equations, one has to replace step by step. For the presented example, the ultimately obtained expressions are

$$\begin{aligned} F &= A^2 \frac{k\_3}{k\_8} \frac{Bk\_5}{Bk\_5 + k\_4} \\ D &= A^2 \frac{k\_3}{Gk\_6} \frac{Bk\_5}{Bk\_5 + k\_4} \\ C &= A^2 \frac{k\_3}{Bk\_5 + k\_4} \\ k\_1 &= Ak\_3r\_1 + \frac{r\_1 + 1}{A} \left( BCk\_5 + BFk\_7 \right) \\ k\_2 &= r\_1 \frac{k\_3A^2 + k\_1A}{B(r\_1 + 1)} \\ k\_0 &= \frac{k\_3A^2 + k\_1A}{r\_1 + 1} \end{aligned}$$

where six parameters are fixed, while one additional parameter r<sup>1</sup> can be chosen freely.

### 3.2. Minus Signs Imply a Low Convergence Rate

For a given data set and a given ODE model, each parameter set determines the time-evolution of the states and its likelihood L can be computed based on the data. Here, parameter values are estimated by minimizing the negative log-likelihood function − log L. For the case of non-linear ODE models, several local optima may exist. In order to find the global optimum, we perform multi-start optimization in combination with a trustregion optimizer. A powerful optimization approach should have a high probability to find a local or the global optimum.

Let us consider an ODE system with four state variables and six reactions of the form

$$\begin{array}{ccc} \emptyset \stackrel{k\_0}{\xrightarrow{k\_0}} A & \qquad \qquad A \stackrel{k\_7}{\xrightarrow{k\_7}} \emptyset & \qquad \qquad A+B \stackrel{k\_2}{\xrightarrow{k\_2}} C\\ \emptyset \stackrel{k\_1}{\xrightarrow{k\_1}} B & \qquad \qquad C \stackrel{k\_3}{\xrightarrow{k\_3}} D & \qquad \qquad D \stackrel{k\_4}{\xrightarrow{k\_4}} \emptyset \end{array}$$

The corresponding ODE system is given by

$$\begin{aligned} \dot{A} &= k\_0 - k\_7 A - k\_2 A B \\ \dot{B} &= k\_1 - k\_2 A B \\ \dot{C} &= k\_2 A B - k\_3 C \\ \dot{D} &= k\_3 C - k\_4 D .\end{aligned}$$

In order to test if negative steady-state expressions lead to optimization problems, we implemented four different steady-state parameterizations, see **Table 5**, and compared the success rate of parameter optimization. For each approach, six parameters are optimized as shown in **Table 5**. Besides the standard approach, i.e., exclusively solving for initial values, two other parameterizations were derived by solving the equation of state B for two different kinetic parameters, namely k<sup>7</sup> and k0. The latter guarantees a non-negative steady-state solution. Apart from that, a fourth parameterization was constructed by adding the equation k<sup>0</sup> = k<sup>1</sup> + 1k<sup>0</sup> to the standard steadystate formulation. In doing so, k<sup>0</sup> is transformed such that k<sup>0</sup> > k1, with the new free parameter 1k<sup>0</sup> describing the difference between k<sup>0</sup> and k1. This approach likewise implies positivity.

For simulation of data, we chose a set of kinetic parameters for ODE integration, initialized the system with its steady-state and excited it by displacement of A at time point t = 30. Data points were generated for 16 different time points by adding normally distributed noise to the model trajectories. In order to study a scenario with different experimental conditions, i.e., different stimulations, simulation was done for three different displacement values, compare cond1, cond2, and cond3 in **Figure 2A**.

**Figure 2A** shows data points and trajectories of a model fit that reached the global optimum. For each steady-state parameterization, we optimized 200 different parameter samples and counted how often several optima were reached. In **Figure 2B**, all converged fits are shown in order of the objective value, in our case the negative log-likelihood value. Several steps corresponding to local optima appear for all steady-state parameterizations, the deepest step corresponds to the global optimum. It can be concluded that the two parameterizations without minus signs, i.e., Solving for k<sup>0</sup> and Standard with Trafo, show a significantly better convergence than the other two. For example, in the parameterization with transformation, the global optimum was twice as often reached than in the standard approach.

In order to explain the convergence behavior of the different steady-state implementations, we analyzed the correlation between initial parameter guess and the success of optimization. The steady-state of the presented model is negative, if and only if k<sup>0</sup> < k1. **Figure 2C** shows the starting samples along the parameter axes of k<sup>0</sup> and k<sup>1</sup> for all four steady-state parameterizations, colors indicate whether a sample did not converge (black) or did converge to a local (blue) or the global optimum (yellow). For comparison, **Figure 2D** shows starting samples along axes of parameters that do not affect the sign of the steady-state, namely k<sup>2</sup> and k4. The sample distribution shows that samples with k<sup>0</sup> > k<sup>1</sup> have a high probability to converge, while samples with k<sup>0</sup> < k<sup>1</sup> tend to abort. On the other hand, the relation of k<sup>2</sup> and k<sup>4</sup> does not have a significant impact on the convergence probability. Furthermore, **Figure 2C** shows that the reparameterized steady-states prohibit sampling in the region with k<sup>0</sup> < k1.

In addition to the starting samples, we analyzed the parameter paths during the optimization. **Figures 2E,F** show the paths for the first 50 starting samples with respect to the above used parameter axes. Parameter samples with k<sup>0</sup> < k<sup>1</sup> usually abort without any considerable steps in parameter space even though several samples cross the border k<sup>0</sup> = k<sup>1</sup> and proceed. In the opposite direction, some samples reach the border when started in the area with k<sup>0</sup> > k<sup>1</sup> and abort exactly at the border. The very most samples drawn with k<sup>0</sup> > k<sup>1</sup> converged to a local or the global optimum G. Again, **Figure 2F** underlines that the convergence behavior is unaffected by the relation of k<sup>2</sup> and k4.

We conclude that steady-state parameterizations that lead to negative parameter values for certain model specifications constitute a severe issue for optimization. Due to the formulation of our steady-state algorithm, negative solutions are automatically avoided in the obtained steady-state expressions.

### 3.3. Dealing with Multiplicity of Steady-States

Let us consider a system with three state variables and seven reactions of the form

$$\begin{array}{ccccc} \emptyset \stackrel{k\_0}{\xrightarrow{k\_0}} A & \qquad \quad A \xrightarrow{k\_1} \emptyset & \qquad \quad \quad \emptyset \xrightarrow{\underbrace{k\_2 \cdot A}{\longrightarrow}} B & \qquad \quad \quad \emptyset \xrightarrow{\xrightarrow{k\_3 \cdot C}} A \\\\ B \xrightarrow{k\_4} \emptyset & \qquad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \xrightarrow{k\_6} \emptyset . \end{array}$$

The production of B is mediated by A, production of A is mediated by C and production of C is mediated by both A and B. The corresponding ODE system is given by

$$\begin{aligned} \dot{A} &= k\_0 + k\_3 C - k\_1 A \\ \dot{B} &= k\_2 A - k\_4 B \\ \dot{C} &= k\_5 A B - k\_6 C .\end{aligned} \tag{18}$$

For this system, the steady-state can still be analytically solved for the states A, B and C. The solution reads

$$\begin{aligned} C\_{1/2} &= \frac{1}{2k\_2k\_3^2k\_5} \left( \frac{1}{k\_1^2 k\_4 k\_6} \Delta + 2k\_0 k\_2 k\_3 k\_5 \pm \sqrt{\Delta} \right) \\ B\_{1/2} &= \frac{1}{2k\_1 k\_3 k\_4 k\_5} \left( k\_1^2 k\_4 k\_6 \pm \sqrt{\Delta} \right) \\ A\_{1/2} &= \frac{1}{2k\_1 k\_2 k\_3 k\_5} \left( k\_1^2 k\_4 k\_6 \pm \sqrt{\Delta} \right) \end{aligned} \tag{19}$$

with the discriminant 1 = k 2 1 k4k<sup>6</sup> · k 2 1 k4k<sup>6</sup> − 4k0k2k3k<sup>5</sup> . For 1 > 0, two positive steady-state solutions S<sup>1</sup> = (A1, B1, C1) and S<sup>2</sup> = (A2, B2, C2) are obtained, while the system has no real steady-state for 1 < 0.

Linear stability analysis reveals that solution S<sup>1</sup> is unstable and solution S<sup>2</sup> is stable, see Supplementary Material. If several steady-state solutions exist, only one of them can be chosen for an optimization run at a time. Here, the stable solution was chosen.

For this system, the issue of multiplicity can easily be solved, however, for more complicated systems several stable solutions might exist. Stable solutions might even switch to unstable solutions along the optimization path, e.g., in case of a Hopf bifurcation. For higher-order equations, analytical solutions become unfeasible and numerical steady-state computation comes into play. However, since the numerical root finding is performed by means of Newton's method, the result depends on initial guesses for A, B and C. Consequently, it is not clear which of the solutions is obtained and stability of the retrieved steady-state is not guaranteed. As we will show, coexistence of stable and unstable steady-state solutions leads to a reduced convergence probability in the numerical approach.

Unlike solving the steady-state equations for A, B and C, the steady-state expressions obtained by our proposed approach are

$$C = A^2 \frac{k\_2 k\_5}{k\_4 k\_6} \qquad B = A \frac{k\_2}{k\_4} \qquad k\_1 = A \frac{k\_2 k\_3 k\_5}{k\_4 k\_6} + \frac{k\_0}{A} \,, \tag{20}$$

where the kinetic parameter k<sup>1</sup> is fixed, while the initial value of A is taken as a free parameter. The obtained solution is unique, since the steady-state equations are linear in the parameters B, C and k1.

In general, our approach avoids multiple-steady-states by choosing a combination of parameters for which the steady-state equations are linear. In doing so, no solution is neglected as long as all steady-state equations are fulfilled. As an analogon, let us consider a single algebraic equation of the form ab<sup>3</sup> + cb<sup>2</sup> + db + f = 0, with the five parameters a, b, c, d and f . On the one hand this equation can be solved for b whereby multiple solutions are obtained. On the other hand it can be solved for one of the parameters a, c, d or f for which the equation is linear leading to a unique solution.

In the following, we compare the convergence behavior of three different steady-state implementations, namely Standard, i.e., analytically solved for the states, Numeric, i.e., numerically solved for the states, and Proposed, i.e., our steady-state approach with positive solutions. For the former two implementations, the seven kinetic parameters k<sup>0</sup> to k<sup>6</sup> are estimated, whereas for the proposed approach, the initial value of A in estimated instead of k1, compare Equation (20).

Since natural systems are always subject to external noise, unstable steady-states are never realized by the system. Therefore, data was simulated by means of the stable steady-state solution. Analogously to Section 3.2, three different displacements of the state A were triggered at time point t = 30 to excite the system. Here, data points were generated for eight different time points.

In order to test the convergence behavior, we started 200 fits from randomly chosen parameter samples. **Figure 3A** shows an example fit that converged to the global optimum. The optimization result of the three steady-state approaches is compared in **Figure 3B**. Steps correspond to local optima. In our approach, nearly half of the samples converged to the global optimum, whereas only about 10% of the fits converged in the standard approach and even less in the numeric approach.