Abstract

Despite of the many experimental studies addressing the evolution and sustainment of heterocyst patterns and the knowledge of the genetic circuit underlying the behavior of single cyanobacterium under nitrogen deprivation, there is still a theoretical gap connecting these two macroscopic and microscopic processes. As an attempt to shed light on this issue, here we explore heterocyst differentiation under the paradigm of systems biology. This framework allows us to formulate the essential dynamical ingredients of the genetic circuit of a single cyanobacterium into a set of differential equations describing the time evolution of the concentrations of the relevant molecular products. As a result, we are able to study the behavior of a single cya-nobacterium under different external conditions, emulating nitrogen deprivation, and simulate the dynamics of cyanobacteria filaments by coupling their respective genetic circuits via molecular diffusion. These two ingredients allow us to understand the principles by which heterocyst patterns can be generated and sustained. In particular, our results point out that, by including both diffusion and noisy external conditions in the computational model, it is possible to reproduce the main features of the formation and sustainment of heterocyst patterns in cyanobacteria filaments as observed experimentally. Finally, we discuss the validity and possible improvements of the model.

Author Summary

When a cyanobacteria filament is deprived of combined nitrogen, some vegetative cells differentiate into heterocysts, Which are terminally differentiated nitrogen-fixing cells. Interestingly, most cells do not differentiate, but remain in their initial vegetative state. The coexistence of heterocysts and vegetative cells is essential for the survival of the filament since (i) heterocysts lose their photosynthetic capacity so they need vegetative cells around to be provided With a source of fixed carbon and (ii) cell division, i.e. reproduction, is only accomplished by vegetative cells. From such a paradigmatic example it is clear that differentiation processes are the result of the interplay of complex regulatory networks acting inside the cell and external stimuli, coming from both the adjacent cells and the environment. In this article we present an integrative approach that combines the study of internal regulatory processes, diffusion, and noisy environments in order to capture the key mechanisms leading to the differentiation of vegetative cyanobacteria into heterocysts and the subsequent pattern formation.

Introduction

The most salient feature that characterizes multicellular organisms is the presence of different cell types, in such a way that the organism associates a different function to each cell type. In each of these cellular types, only a subset of the genes that constitute the genome of the organism (genotype) are expressed, which identify the function and morphology of the cell (phenotype). The development of specialized cells involves differentiation processes, which lead to alterations in gene expression producing different phenotypes from a given genotype. These processes are highly dynamical, directed by complex regulatory networks involving cell-to-cell interactions, and often triggered by external stimuli. As a result of the differentiation processes a rich cooperative pattern involving different cell types is established, increasing the complexity and adaptability of the organism. Due to the large number of scales involved, ranging from protein binding to diffusion of specific elements throughout the organism, a correct mathematical modeling of differentiation processes and their associated pattern formation demands an integrative approach combining tools from statistical mechanics and the theory of dynamical systems (see [1, 2] for instance).

Cyanobacteria are one of the first organisms that developed multicellularity some (2—3) billion years ago [5]. These bacteria perform oxygenic photosynthesis releasing oxygen to the environment. However, nitroge-nase, the enzyme that performs nitrogen fixation, is deactivated by oxygen so that nitrogen fixation cannot occur in its presence [6]. Cyanobacteria solve the incompatibility of incorporating both oxygenic photosynthesis and nitrogen fixation by separating these processes (1') temporally, such as in the unicellular Cyanothece sp. strain ATCC 51142, which presents photosynthetic activity during the day and fixes nitrogen during the night [7] , or (ii) spatially, by the generation of non-photosynthetic nitrogen-f1xing cells distributed along the filament and acting as nitrogen suppliers.

However, in the absence of combined nitrogen (cN), a subset of the vegetative cells differentiate into heterocysts, which are terminally differentiated nitrogen-fixing cells. By differentiating, heterocysts lose their photosynthetic capacity, so they require an external source of fixed carbon [8, 9]. To this aim, each forming heterocyst sends a signal, by means of the production of some substance that diffuses along the filament, to prevent the differentiation of its neighboring cells. A cooperative pattern is thus established: heterocysts provide cN to the filament while vegetative cells supply fixed carbon. As a result, heterocysts appear interspersed with around 10 vegetative cells, depending on the species, forming a semi-regular pattern that remains approximately constant regardless of cell division [10, 11]. The resulting pattern forms one of the simplest and most primitive examples of a multicellular organism as a product of the interdependence between heterocysts and vegetative cells. Interestingly, an isolated cyanobacterium does not differentiate but it first divides so that one of the descendants differentiates. This latter mechanism is crucial since (i) a sole heterocyst would lack a source of fixed carbon and (ii) it would not reproduce as it is a terminally differentiated cell [12].

In references [13, 14] Rutenberg and coworkers analyzed a model to explain heterocyst patterns by means of the study of cN diffusion along a cyanobacterial filament. On the other hand, Gerdtzen et al. [15] modeled cyanobacterial filaments based on a time-discrete dynamical system incorporating the main interactions between the most important proteins that take part in heterocyst formation. In this work, we develop a simple mathematical model by incorporating the recent experimental results on the genetic regulatory network of cyanobacteria into the theoretical machinery of system biology.

Furthermore, our model shows that noise plays an important role in the onset of differentiation by enabling the development of the characteristic heterocyst patterns for a wide range of model parameters. This reveals that cyanobacteria filaments have developed an efficient response to the noisy conditions that characterize the natural environment.

First we present the main actors of the basic regulatory network and the different dynamical interactions that take place during the differentiation process. Then we develop a mathematical model for the unicellular reaction to nitrogen deprivation. Although a single cell model cannot provide a complete understanding of heterocyst formation, we analyze the main features that arise from the dynamical behavior of the system to gain insight about cell dynamics under different external conditions. Finally, we round off the paper by introducing the spatial model consisting of a filament of cyanobacteria, each one characterized by the dynamical circuit developed previously, that interact by means of protein diffusion.

Results

Description of the main genes and their basic genetic circuit

This process is usually completed after 20 hours at 30°C [9]. In Fig. 1 we show a basic scheme of the genetic circuit including the most relevant elements and their respective interactions. Here we explain the main features of this genetic circuit. more, the transcription of the genes targeted by NtcA does not start in the absence of 2-OG [22, 23].

CYANOBACTERIUM

Main components and interactions involved in the reaction to combined nitrogen deprivation in cyanobacteria. Rectangular boxes represent genes (ntcA, bet]? and patS) while rounded boxes and circles represent transcription factors (NtcA, HetR and PatS) and smaller molecules (2-OG and 0N) respectively. Normal-tipped and flat-tipped arrows stand for up-regulating and down-regulating processes respectively. Dashed lines stand for indirect or imperfectly understood interactions. The accumulation of 2—OG enhances the DNA-binding activity of NtcA, which in turn up-regulates the transcription of ntcA and hetR. HetR activates ntcA and bet]? (composing the central NtcA-HetR autoregulatory loop), the inhibitor pats and other genes that lead to nitrogen fixation and the morphological changes involved in heterocyst differentiation. 2—OG and 0N levels are linked through the GS/GOGAT cycle (see Fig. 2).

To bind DNA, NtcA needs to homodi-merize [29, 30]. In conclusion, the accumulation of 2-OG is the factor that triggers differentiation. In agreement With this idea, artificial increased levels of 2-OG result in heterocyst development even in the presence of ammonium [16, 18, 31].

Remarkably, null mutants of hetR do not produce heterocysts Whereas an overeXpression of hetR leads to an increased heterocyst frequency [27, 32, 33]. The transcription of hetR is induced by NtcA through the action of an intermediate, nrrA [28]. The DNA-binding activity of HetR requires its homodimer-ization [34, 35]. Multiple transcription factors related to heterocyst formation are up-regulated by HetR, including hetR itself [35], ntcA [36] and patS [35].

However, the action of NtcA and HetR alone cannot eXplain pattern formation. Another transcription factor, PatS, inhibits the DNA-binding activity of HetR [8, 35, 40, 41]. This inhibitory behavior is essential for the communication With adjacent cells and thus to achieve the observed patterns of vegetative cells and heterocysts in cyanobacteria filaments (see Fig. 3). Furthermore, patS is strongly eXpressed in differentiating cells and mature heterocysts due to its 2-Oxogluta rate

GOGAT

Glutamate Glutamine

GS

GSIGOGAT cycle. 2—OG and cN indirectly interact through the GS/GOGAT cycle. Glutamine is transformed into glutamate by means of 2—OG through the 2-OG amidotransferase (GOGAT) while cN converts glutamate into glutamine through the glutamine synthetase (G8). The importance of the cycle in heterocyst differentiation is twofold. From one side, it constitutes the early one-cell sensor to nitrogen starvation: the absence of cN breaks the cycle down and 2—OG starts to accumulate, whose action leads to the cascade of processes that provoke the differentiation (see Fig. 1). Additionally, later during the differentiation, it processes the cN created by the heterocysts decreasing the levels of 2—OG. The latter is crucial forthe formation of the heterocyst pattern (see Fig. 3).

A filament Without patS develops multiple contiguous heterocysts (about a 30% of all cells as compared to the usual 10% in the Wild-type filament). On the other hand, an over-eXpression of patS suppresses heterocyst differentiation [9]. Moreover, the addition to the growth medium of a synthetic peptide composed of the last five residues (RGSGR) of PatS (PatS5) inhibits heterocyst development, suggesting that PatS5 may be a diffusive mature form of PatS that stops the differentiation of the rest of vegetative cells of the filament

m Vegetative

Diffusion scheme. Schematic representation of the diffusion processes that sustain the heterocyst pattern. Heterocysts produce cN and PatS. cN diffuses along the filament where, due to the action of the GS/ GOGAT cycle (see Fig. 2), decreases the levels of 2—OG breaking the autoregulatory core NtcA-HetR. Early during the differentiation, PatS (or other derivative of it, see the text) diffuses along the filament inhibiting HetR. Both processes combined prevent the differentiation of the rest of vegetative cells and explain the formation of the pattern.

A similar protein carrying the RGSGR pentapeptide, with a similar effect as that of PatS, is HetN [42]. A chain lacking both PatS and HetN leads to a lethal phenotype in which all cells differentiate [43].

To this end, two new membrane layers are biosynthesized to decrease the entry of oxygen into the cell [44]. The morphogenesis of these two layers is controlled by two family of genes, hep and hgl, that are indirectly up-regulated by HetR [35]. After these morphological changes the genes in charge of nitrogen fixation, genes, are expressed. These genes encode, among others, the enzyme ni-trogenase, which ultimately performs nitrogen fixation.

Thus, the diffusion of these inhibitors from heterocysts along the filament plays a key role in pattern maintenance (see Fig. 3). As a result of the differentiation het-erocysts produce f1xed nitrogen from N2 of the atmosphere and they interchange this nitrogen with the oxygen derivatives produced by the vegetative cells, in an illustration of the cooperative behavior between cell types in a multicellular organism.

Regulatory equations

Details are left to supplementary information (Sl text). To simplify notation, constants related to NtcA, HetR, PatS and cN are denoted With the letters a, r, s, and 11 respectively.

We assume that the probability that NtcA binds the promoter in the absence of 2-OG can be neglected. Taking into account that both HetR and NtcA dimerize to bind DNA we find:

La, the so-called leak term, measures the basal production of ntcA in the absence of regulation. Subscripts and superscripts identify the binding site and the transcription factor for Which the constants are given respectively.

We assume that hetR is regulated by NtcA by means of a usual Hill function, yet the real process presents an intermediate, nrrA. To do so, we take into account that nrrA concentration relaxes rapidly to a limiting value. Furthermore, PatS affects the autoregulatory loop of HetR. It has been suggested that PatS binds the binding site of HetR in the promoter of hetR preventing HetR binding [35]. These facts, along With the influence of 2-OG levels, provide With an expression for the transcription velocity: HetR regulates most processes of the genetic circuit. It governs, among others, the transcription of ntcA, patS, hep, hgl and genes that lead to most of structural changes of the cell and to nitrogen fixation.

For simplicity, we implicitly include the effect of HetN in PatS, as their action is expected to be equivalent (see the previous section). This gives the simple transcription velocity:

Let us begin by examining nitrogenase concentration, which is directly controlled by genes. Although this is not a direct process, we can assume, as we did for the NtcA-regulation of hetR, that genes are functionally governed by [HetR] following a typical Hill function. The nitrogenase production rate is given by:

We can effectively account for the lag introduced by intermediate processes not taken into account explicitly in the model by increasing the value of 6N1 so that [Ni] relaxes more slowly. Assuming that nitrogenase produces fixed nitrogen at a constant rate, we arrive at the equation that governs cN levels in cyanobacteria:

Assuming that the levels of cN relax rapidly we solve Eq. (5) for the steady state. Substituting in Eq. (4) we find: Where

Both are related by means of the GS/GOGAT cycle (Fig. 2). Assuming the cycle is in equilibrium and reactions are grounded on the law of mass action, the following two conditions must be satisfied: Which lead to the relation:

In fact, 2-OG production is controlled by some processes that are not considered in this work and so its value must be limited. We can effectively include such a limiting value by means of a translation on

They represent the temporal variation of the most important factors of the genetic circuit, namely NtcA, HetR, PatS and cN. Using the production rates (1), (2), (3), (6) and introducing degradation rates constants, 6*, we find: Where we have introduced the dimensionless variables: and the constants

The study of the cyanobacterial filament is left to the final section. We show that the main modification Will be adding diffusion processes for the inhibitors PatS and cN through the chain. An important ingredient in pattern formation, noise, Will be also added to the equations.

Unicellular dynamics

In this section we analyze the dynamical system (11) for a set of constants (Table 1) that exhibit both the dynamical and the structural properties of heterocyst differentiation. Following the usual practice in the analysis of dynamical systems, we study the basic properties of equations (11), such as fixed points and linear stability analysis, to analyze the key features leading to heterocyst differentiation. Taking into account the difference between the relaxation times of the constituents of the model, given by the inverses of d* (see Table 1), we can interpret it as composed of two temporally separated systems: a rapid one, formed by HetR and NtcA, showing fast dynamics that relaxes to its steady state almost instantaneously and a slow one, composed of PatS and cN, whose evolution is dictated by the values of HetR and thA in their instantaneous equilibrium.

Parameters for Eq. (11) that reproduce heterocyst formation under noisy conditions and pattern formation when Fats and cN diffuse along a filament of cyanobacteria. This corresponds to an adiabatic elimination technique [49] that helps in reducing the complexity of the dynamical system by splitting it into two simpler interdependent subsystems. First, we look at the fixed points of qa and qr for each pair of values of qs and q”

We find three different branches of solutions that coexist in some regions. The fixed points on the lower and upper branches are always stable (blue region in Fig. 41) and those lying on the middle branch (red region) are saddles. Transitions between the regions with one and three fixed points correspond to saddle-node bifurcations in which the middle branch of solutions coalesce with the lower and the upper one respectively. The basins of attraction of both stable fixed points are separated by the stable manifold of the saddle point (Fig. 411).

A sufficient large perturbation may result in the system crossing the manifold of the saddle and falling into the other stable branch of solutions. The distance between the saddle and the nodes determines the size of the perturbation needed to activate or inactivate the system.

In the regions showing bistability the field takes two very different forms, one corresponding to the values of the lower branch and another corresponding to those of the upper one (Fig. 5). We expect a hysteresis effect: if initially the dynamics lies on a particular branch it will remain on it unless a fluctuation or a bifurcation makes the system jump to the other branch.

5A) we find only one stable fixed point that corresponds to a vegetative state (lower branch). The upper branch is completely unstable: any dynamics lying on it will fall down to the lower branch and eventually be attracted to the vegetative sink. The steady state is very robust against perturbations since it is far from the bifurcation region and there is a significant distance to the saddle in the qr — qa plane. By reducing the flow of cN from the exterior of the cell (ln = 0) we find that a stable fixed point appears in the upper branch, a heterocyst state, while the vegetative state gets closer to the bifurcation region, thus becoming more susceptible to perturbations that can make the system reach the upper branch. In the absence of cN, the cyanobacterium would evolve from state A to state B in the lower branch until a perturbation pushes it to the upper branch, eventually becoming an heterocyst due to the field acting on that branch (Fig. 5B and C).

The heterocyst fixed point also becomes more stable due to diffusion, since its production of inhibitors is distributed among other cells (see Fig. 5D).

Strains of cyanobacteria.

Heterocyst patterns

There we have shown that, for a specific range of parameters, the model exhibits features that would lead to heterocyst development under noisy conditions. Nevertheless, the model should be extended to cyanobacteria chains to account for heterocyst development since, as previously noted, isolated cyanobacteria do not become heterocysts by themselves; the action of the chain is needed to generate heterocysts.

The main modification is the introduction of diffusion of PatS and cN along the cyanobacteria chain. For this purpose we add to Eq. (12) the discrete version of the diffusion equation: Where DC is called the difiusion constant of the element C. Now, it is straightforward to introduce PatS and cN diffusion into the equations. The dynamics of cell i is characterized by the following set of equations: which constitutes the model for a cyanobacteria filament. To account for environment variability we add white noise, G; *(t), of the same amplitude, (G; *(t) G; *(tj)) = 560‘ — 1’), for all the components of the system. Based on these equations, we investigate the conditions that lead to a heterocyst pattern. It is easy to notice that they correspond to an activator-inhibitor system of cells coupled in a reaction-diffusion scheme [50]. This kind of system produces regular pattern formation [51—53]. Turing (linear stability) analysis of equations (16) (see 82 text) provides insight on the periodicity of patterns. It is interesting to show that the minimum periodicity observed in such analysis is larger than 1, which means that a single bacteria is unable to differentiate. We performed the direct integration of equations (16) for chains of 200 cyanobacteria. We used a Runge-Kutta method for the numerical integration of stochastic differential equations

All simulations were performed with periodic boundary conditions, i.e. emulating a circular filament, for simplicity. We have also tested the more realistic no flux boundary conditions and find no change in gene dynamics and heterocyst patterns in the interior of the filament. This shows that the effect of boundary conditions is highly localized around the borders. In simulations with no flux boundary conditions no heterocysts were found in the border, which is in good agreement with experimental observations [55]. The level of noise that best reproduces heterocyst pattern is f = 0.001 for the set of parameters of Table 1. Importantly, isolated cells do not initiate differentiation with this level of noise, in agreement with the results from the linear stability analysis. Diffusion constants have been set to D5 = 0.1 and D” = 0.2. Heterocysts patterns develop for different levels of noise and diffusion constants, but the model parameters, which characterize cell response to nitrogen deprivation, should change accordingly. This correlation between noise, diffusion and model parameters supports the idea that cyanobacteria have evolved towards a better response to the normal levels of noise in their environment.

We observe that the filament concentrations relax to the constant protein levels of the vegetative state we showed in the previous chapter. Then, due to the coupled action of noise and diffusion, some cells start to differentiate. As new forming heterocysts appear, their production and exportation of inhibitors to the surrounding cells make the latter more stable to perturbations stopping their differentiation. The model reproduces very well the initial peak that both NtcA and HetR present experimentally [35, 56]. PatS increases more slowly to its steady value reducing the levels of NtcA and HetR and, finally, cN is generated by heterocysts stabilizing the pattern.

7 shows the evolution of the profile for a 200 cells chain of cyanobacteria. We observe that heterocysts progressively appear in those regions in which other heterocysts do not have effect (1'. e. those vegetative cells that are not supplied of sufficient cN and PatS). Finally, a semi-regular pattern is generated. PatS and cN diffuse along the filament exhibiting smooth variations between vegetative cells and heterocysts, while HetR and NtcA present very abrupt variations between cell types.

It should be stressed that although initially some close heterocysts appear, they are eliminated by the nonlinear action of the system during the differentiation process. Close heterocysts compete for the same region of action (the same vegetative cells that consume their PatS and cN) and then they cannot reach the optimal heterocyst state, which is stable to slight perturbations. Finally, one of them falls down from the upper branch becoming a vegetative cell. This behavior is typically observed experimentally [4, 9]. The final histogram can be nicely fitted by a F-distribution, implying Poisson-distributed wait times for the underlying noise driven process.

Discussion

This phenomenon is the basis for multicellular organism and pattern formation. The approach presented here deals With a simple system, heterocyst formation in cyanobacteria filaments, yet compleX enough to capture the main ingredients of some of the mechanisms for cell differentiation and pattern formation under external driving. The knowledge of the basic regulatory genes and their corresponding interactions allows for a detailed description of cell dynamics. We have derived the evolution equations of the involved genes based on the statistical mechanics of their corresponding regulatory processes. This allows us to obtain a detailed description of the continuous time dynamics of the main regulatory proteins, in contrast to other discrete approximations based in boolean dynamics [15]. This kind of analysis has been shown to describe successfully other time dependent phenomena concerning cyanobacteria such as their circadian cycles [57]. The analysis of the unicellular dynamics has revealed that the two cellular stable states, vegetative and heterocyst cells, appear as attrac-tors of the nonlinear dynamics of the regulatory equations. However, the study of many coupled cells is needed as cyanobacteria do not differentiate when isolated.

We have shown that one important ingredient affecting the dynamical behavior of the chain is noise, which plays a key role in onset of the pattern formation, i.e., the transition from the initial chain of vegetative cells to the steady state in which heterocysts coexist with vegetative cyanobacteria. Thus, the appearance of differentiation is, in our model, a pure stochastic event. The cooperative character of the filament is clear from the amount of noise needed to start the differentiation process which appears significantly smaller than that needed in isolated cells. The source of noise as well as its biological consequences is, nowadays a current topic of research [58]. In fact, at its initial state, differentiation of cells appears randomly along the filament, but shortly after its onset a characteristic distribution of heterocyst emerges. This distribution can be compared with the experimental one with a fairly good agreement [41].

One issue that have not been considered in this work is the replication of vegetative cells. This effect has been taken into account in [14]. Although this improvement is relevant, it only affects, in our approach, to the mean separation between heterocysts, by opening a gap in the F-function shape of Fig. 7B and thus approaching better to the experimental distribution.

Unlike other approaches [13] in which comparison is done (globally) with heterocyst distributions, our work would allow for a qualitative comparison of each component involved in the differentiation (see Fig. 6). Unfortunately, there is not enough experimental data to make a detailed fit so to extract reliable parameters. The availability of such data is extremely important both for having a better set of model parameters and to validate new models. A complete understanding of the mechanism that derive in phenotypic differentiation is the first step for a modular comprehension of the whole cell [59].

Methods

Eq. (16) is a set of stochastic differential equations (SDE) so its numerical integration requires generating a statistical representative trajectory for a discrete set of time-values. A SDE of the form Where G(t) is a Gaussian White noise With can be integrated through a Runge-Kutta integration algorithm by adding a particular Gaussian signal at each stage of the scheme. This algorithm coincides With the usual Runge-Kutta scheme for f = O. In this work we have employed a 30432G algorithm, Which is correct up to 3th order, is developed in 4 stages and uses 2 independent Gaussian random variables.

Supporting Information

Regulatory equations: a statistical mechanics approach. (PDF)

Turing linear stability analysis. (PDF)

Acknowledgments

We are grateful to S. Ares and I. Munoz-Garcia for sharing ideas and useful discussions. We also acknowledge the “Genetic Regulation and Physiology of Cyanobacteria” group at University of Zaragoza for sharing insight on the genetic of heterocyst formation.

Author Contributions

Performed the experiments: ATS. Analyzed the data: ATS IGG FF. Contributed reagents/materials/analysis tools: ATS IGG FF. Wrote the paper: ATS IGG FF.