Modeling heterocyst pattern formation in cyanobacteria
© Gerdtzen et al. 2009
Published: 16 June 2009
Skip to main content
© Gerdtzen et al. 2009
Published: 16 June 2009
To allow the survival of the population in the absence of nitrogen, some cyanobacteria strains have developed the capability of differentiating into nitrogen fixing cells, forming a characteristic pattern. In this paper, the process by which cyanobacteria differentiates from vegetative cells into heterocysts in the absence of nitrogen and the elements of the gene network involved that allow the formation of such a pattern are investigated.
A simple gene network model, which represents the complexity of the differentiation process, and the role of all variables involved in this cellular process is proposed. Specific characteristics and details of the system's behavior such as transcript profiles for ntcA, hetR and patS between consecutive heterocysts were studied.
The proposed model is able to capture one of the most distinctive features of this system: a characteristic distance of 10 cells between two heterocysts, with a small standard deviation according to experimental variability. The system's response to knock-out and over-expression of patS and hetR was simulated in order to validate the proposed model against experimental observations. In all cases, simulations show good agreement with reported experimental results.
A simple evolution mathematical model based on the gene network involved in heterocyst differentiation was proposed. The behavior of the biological system naturally emerges from the network and the model is able to capture the spacing pattern observed in heterocyst differentiation, as well as the effect of external perturbations such as nitrogen deprivation, gene knock-out and over-expression without specific parameter fitting.
Cyanobacteria are blue-green algae, prokaryotic organisms capable of both obtaining energy by oxygenic photosynthesis and nitrogen fixation. They play a fundamental role in earth's carbon cycle as primary producers, and in the nitrogen cycle as nitrogen fixers . These two processes are somewhat contradictory in nature, since the activity of the nitrogen fixing enzyme, nitrogenase, is dramatically reduced by the presence of oxygen, which is the end product of photosynthesis. Cyanobacteria evolved a mechanism for fixing dinitrogen (N2), an oxygen sensitive process, under aerobic conditions. Vegetative cells, which are the cells where oxygen is produced from photosynthesis, can differentiate into specialized non-dividing nitrogen-fixing cells called heterocysts. Anabaena and other filamentous cyanobacteria in the order of Nostocales grow in filamentous structures formed by photosynthetic vegetative cells. For some of these strains, under nitrogen-limiting conditions, vegetative cells differentiate into heterocysts at intervals along the filaments therefore generating a semi-regular spacing pattern . This mechanism has made Anabaena a model organism for studies involving cellular differentiation and pattern formation in prokaryotes.
The primary known function of heterocysts is the fixation of dinitrogen, for which they require a reductant supplied by the vegetative cells since differentiated cells gain the ability to reduce nitrogen gas but lose the ability to fix carbon dioxide. Heterocysts have a deactivated photosystem II (a protective envelope that is semi permeable to oxygen) and higher respiration rates in order to consume the oxygen that was able to enter the cell. All these characteristics generate an anaerobic environment suitable for the operation of the nitrogenase enzyme and therefore, for nitrogen fixation [3, 4]. Nitrogen fixed in heterocysts is transported along the filament, likely in the form of glutamine, and utilized by the whole population of cells . In this scheme, heterocysts share nitrogen with vegetative cells receiving in turn carbon resources establishing a cooperative system .
Differentiation of vegetative cells into heterocysts is triggered by the absence of a fixed nitrogen source in the growth medium . A large number of genes involved in the development and spacing of heterocysts have been identified, however the complete mechanism by which they interact is still unclear . Interaction mechanisms, based on current knowledge and inferred interactions have been proposed [8–10].
The reduction in nitrogen levels rapidly enhances the activation of ntcA, a DNA-binding factor involved in the transcription of genes involved in nitrate and ammonium transport and assimilation, dinitrogen fixation and heterocyst development, including hetR, which is required for heterocyst differentiation [11, 12]. Nitrogen levels are sensed by the cell's intracellular levels of 2-oxoglutarate. This intermediate from the Krebs Cycle accumulates in the cytoplasm upon nitrogen starvation, enhancing the DNA binding activity of NtcA . There is evidence that NtcA binds to the promoter region of its own gene, suggesting that it regulates its own expression .
hetR plays a key role in the regulation of heterocyst differentiation. This gene appears to be indirectly activated by NtcA through intermediate signaling molecules, since none of its promoters contain the sequence to which NtcA binds . The gene product HetR is a DNA binding protein that acts as a homodimer . In the presence of nitrogen, hetR is transcribed at low levels in all cells. Following deprivation, hetR expression is induced in a process that requires ntcA expression and the presence of a functional hetR gene, indicating that hetR is also autoregulatory .
For the generation of the characteristic pattern observed upon heterocyst differentiation, besides the activation of genes required for heterocyst development, the genes responsible for regulating the position and spacing of heterocysts must be activated. Among the genes involved in the control of heterocyst spacing are hetN, patA, patS and a few others. The patA gene is transcribed at low levels in the presence of nitrogen, and transcription increases following nitrogen step-down in a similar pattern of expression to that of hetR. PatA is a protein whose structure suggests kinase activity, which is speculated to act indirectly on controlling heterocyst spacing since it does not contain a known DNA-binding motif. This protein is required for the enhancement of hetR transcription in developing heterocysts [9, 15]. HetN regulates the activity of hetR likely through a secondary metabolite, therefore playing a regulatory role in the control of heterocyst formation .
PatS acts as a spacing regulator for heterocyst differentiation. This protein has been shown to block heterocyst differentiation when present in multiple copies or when over-expressed. PatS is produced early in heterocyst development and is believed to be processed (maybe by HetR) to release a pentapeptide capable of diffusing through the filament's periplasmic space to form an inhibition gradient around the developing heterocysts. This is supported by the finding that a synthetic peptide corresponding to the last five amino acids of PatS inhibits heterocyst development . Studies have shown the feasibility of diffusion of this pentapeptide through nonspecific intercellular channels from cytoplasm to cytoplasm in Anabaena [20, 21]. In addition, when patS is knocked-out filaments show an increased frequency of heterocysts and an abnormal spacing pattern . It is speculated that PatS producing cells must be refractory to the influence of PatS, in order not to affect hetR gene expression, since it seems paradoxical how HetR can promote its own synthesis in the presence of the inhibitory pentapeptide derived from PatS. In fact it has been proposed that PatS is inactive in the cytoplasm of heterocysts and that an active peptide is released into the periplasm and transported to vegetative cells where it would act on HetR, or act on the transcription factors directly involved in the expression of hetR .
Several mathematical models have been developed in order to systematically explore and understand the system's behavior. Baker and Herman proposed a mathematical model for differentiation where heterocysts produce a substance that diffuses along the filament and inhibits the differentiation of vegetative cells into heterocysts. When its level in a cell drops below a threshold value, the cell differentiates into a heterocyst. The observed pattern of heterocysts is reproduced but the model shows high sensitivity to the parameters chosen, many of which are vague in terms of association to biological phenomena . This model involves the production of an inhibitor by heterocysts and its diffusion along the filament where heterocysts develop if the inhibitors' concentration drops below a certain threshold. A subsequent work explains heterocyst differentiation by means of a competitive mechanism [23, 24]. This mechanism involves an autocatalytic activator of differentiation produced by proheterocysts and an inhibitor produced in proportion to their degree of differentiation, which stops and even regresses proheterocyst differentiation and whose effect is reduced by adjacent vegetative cells. With this mechanism an irregular pattern will form on a filament where a maximum and minimum distance between heterocysts is observed . Some of the elements in this proposed mechanism are in agreement with the underlying gene network currently identified to control differentiation of vegetative cells into heterocysts and addressed in this work. In this network the differentiation activator can be associated with the protein HetR, which enhances its own production, and the inhibitor with PatS, which inhibits the production of HetR .
Turing developed models that showed the possibility of obtaining patterns from an initially homogeneous distribution, from a reaction-diffusion system where two substances with different diffusion rates interact . Gierer and Meinhardt showed that in morphogenesis a spontaneous pattern formation is possible if a locally short-range self-enhancing reaction associated with an activator which promotes its own production, is coupled with an antagonist that acts on a longer range [27, 28]. Using this hypothesis Meinhardt showed that the Anabaena system satisfies the prediction of this model, considering a locally restricted activator HetR that self-activates nonlinearly (as a dimer) and a diffusible inhibitor PatS . The pattern formed by this network has also been modeled and captured using L-system models where the network is described as a growing system with reaction-diffusion equations for antagonic activator and inhibitor molecules . A population model was proposed by Pinzon and Ju to describe the effects of cellular activities and cultivation conditions on heterocyst differentiation at a culture-level . More recently, Allard et al. proposed a dynamic model that considers random cell growth and division, as well as production, transport, and consumption of fixed nitrogen within the filament. In this model, nitrogen is the main trigger of cell differentiation and cell growth acts as a stochastic factor that can induce cell differentiation. Simulation results are in good agreement with experimental pattern distributions .
In this paper, a simple evolution mathematical model is presented, to study the pattern formation obtained upon heterocyst differentiation, based on the main regulatory elements and interactions of the gene circuit involved in differentiation: ntcA, hetR and patS. In contrast, previous models of the system are based on the protein network involved in heterocyst differentiation. The mathematical model proposed behaves as a discrete one-dimensional activator – inhibitor system, analogous to a continuous autocatalysis-inhibition model [26, 27], but directly derived from the genetic network of the cyanobacteria. This model allows, not only to capture the stable states obtained by the system upon differentiation, but also the response of the system following modifications such as knock-out or over-expression of some of the elements and their effect on the patterns formed as well as the stable states achieved by the system. The underlying idea is the same as in previous models: there is an emerging periodic structure of localized heterocysts in the presence of a local autocatalysis (HetR) coupled with and a long-range inhibition (PatS). Nevertheless, the discrete model presented here has several important advantages when compared with the classical continuous autocatalysis-inhibition models. On one hand, the proposed model naturally arises from the interactions of the gene network involved in differentiation with a set of parameters that have direct biological interpretation. This is a bottom-up construction alternative to the typical top-down methodology of the continuous case where the parameters have to be adapted to each case from general set-up models.
Once the discrete model is established, it allows computing the patterns of the system upon differentiation and the frequency of the corresponding attractor basins. These are the equilibrium states of the system reached from a representative sample of random initial conditions. These types of calculations are a difficult task for the system of partial differential equations in continuous models. Moreover, since the model parameters have a biological meaning, the effect of modifications on the genetic network, such as knock-out or over-expression on the patterns achieved, can be tested. For this reason, the proposed model is suitable for experimental design and for exploring the pattern formation capabilities of systems that could have these pattern forming mechanisms. Finally, some variants of the model that include cell growth and death, or reduced binary on-off activity states can be implemented, but were not considered here for the sake of simplicity.
A finite array of n genetically identical cells has been considered. A fixed number of n cells in the array, organized in a cyclical manner where each cell is connected to its right and left neighbors is used.
where the values of row i (i = ntcA, hetR, patS) represent the activation or inhibition strength that the proteins of the intracellular gene circuit exert on i. For instance, no protein acts on ntcA intracellularly (row 1) and ntcA acts on hetR with an activating strength of one. C Int summarizes all the network interactions that occur intracellularly for a single cell.
A reduction in nitrogen availability triggers the initiation of heterocyst differentiation mechanisms. Initially all cells behave in the same way, and the internal activity of each individual cell causes an inhibitory effect on two elements of the intracellular network of cells in its immediate neighborhood. Low external nitrogen levels trigger an increase in ntcA expression. NtcA leads to the transcription of hetR, which is self-activating. The increase in HetR levels enhances patS transcription to a lower extent than the other activation processes occurring intracellularly (represented with a strength of 0.5 in C Int ), resulting in higher levels of the process pentapeptide PatS-5 (represented by PatS in Figure 1). It has been shown that this pentapeptide has an indirect inhibitory feedback effect on the expression of hetR in adjacent cells, by preventing the DNA binding of HetR and subsequent hetR up-regulation, therefore blocking heterocyst formation [16, 19]. This is represented by a -1 in C Int . There is also an indirect long-term effect produced by the increased amount of nitrogen that becomes available once vegetative cells differentiate into heterocysts and start fixing nitrogen. This (indirect) effect is induced by HetR and sensed by the cell as low oxoglutarate levels, which in turn decrease the binding activity of NtcA . Both PatS and the products of nitrogen fixation are the main signals determining the heterocyst pattern formed .
This matrix shows the fact that ntcA expression is reduced by the indirect extracellular action of HetR and hetR expression is reduced by the action of extracellular PatS. C Ext summarizes all the interactions that occur in the network for a single cell due to its extracellular environment.
resulting a circulant (Toeplitz) type matrix A. The global interaction matrix A summarizes all the interactions that occur both intracellularly and extracellularly for an array of n cells with a gene network as the one shown in Figure 1.
The state of the system is associated with the protein levels of each of the elements in the network (ntcA, hetR and patS). These protein levels are normalized and associated to a real number between 0 and 1. The state of the system is therefore described by the levels of these factors on each cell of the array, represented by the vector x∈ ℝ3n .
with u 0 = 1 in the absence of nitrogen and u 0 = 0 in its presence. Each triplet (1 0 0) of the vector uis associated to the effect of an external input on an element of the network (ntcA, hetR, patS) for a cell in the array. A value of 1 indicates that nitrogen levels directly affect the first element on the network for each cell, ntcA.
The absolute value of deviation of from 10, , for several values of D is shown in Figure 3-B. A lower value indicates a closer to 10. In this case 5000 biologically feasible random initial conditions were used and D was sampled in the interval [0.730; 0.810] every 0.001 units. It was found that the transport factor (D 10) that minimizes the deviation DEV 10 in Figure 3-B is 0.767 for systems with 50 to 100 cells. Using this D 10 gives systems with 's only 0.12% far from 10 (average error). The average of those transport factors is = 0.767 ± 0.001. Based on this fact, 0.767 will be used as the transport factor for the following simulations.
Average and standard deviation for the distance between heterocysts in the absence of nitrogen for the wild type , when patS expression has been knocked-out and when hetR has been over-expressed . Simulations consider a transport factor D = = 0.767 and 5000 biologically feasible randomly chosen initial conditions.
Number of cells
10.010 ± 2.357
4.603 ± 3.368
6.161 ± 3.730
9.991 ± 2.356
4.611 ± 3.382
6.135 ± 3.725
9.977 ± 2.353
4.627 ± 3.382
6.119 ± 3.733
10.007 ± 2.355
4.627 ± 3.397
6.164 ± 3.727
9.999 ± 2.364
4.627 ± 3.398
6.115 ± 3.721
9.979 ± 2.364
4.612 ± 3.386
6.084 ± 3.691
The maximum ntcA expression levels change with L H and cell array size. In fact, as the distance between heterocysts increases, ntcA expression increases, as corroborated by experimental observations. Similarly, the inhibitory action of neighboring cells on ntcA expression is reduced as the distance between heterocysts increases. If the distance between heterocysts is high enough, oxoglutarate levels will be high to enhance ntcA expression on the cell located at position L H /2. Eventually, this expression will be high enough to enhance hetR expression and trigger cell differentiation into a heterocyst.
The differentiating behavior of vegetative cells into heterocysts is triggered by external disturbances on nitrogen levels, and the cellular response to these disturbances based on its internal characteristics. In the presence of nitrogen, ntcA expression is reduced and differentiation does not occur. On the other hand, nitrogen absence will trigger ntcA expression and the chain of interactions leading to the differentiation of some cells (see Figure 1) .
The model presented in this paper is able to capture both of these behaviors, as has been shown in the previous sections. As an illustration, a case study is shown.
It has been reported that genetic manipulation of the gene network involved in cell differentiation into heterocysts does change the differentiation pattern observed. Two extreme examples are the knock-out and over-expression of patS . Simulation of the differentiation of an array of 100–50 cells with D = with knocked-out expression of patS was performed. 5000 randomly chosen initial conditions were considered. Figure 6-C shows that for patS knock-outs the inhibitory effect of PatS is reduced, therefore decreasing the average distance between heterocysts and even allowing the formation of heterocyst clusters in some cases. Results in Figure 6-C are supported by L H histograms shown in Figure 4-B. These results indicate that knocking-out the expression of patS shifts L H distribution, significantly shortening the average distance between heterocysts, , compared to the wild type (see Figure 4-A), as observed experimentally . This is consistent with the results presented in a previous model, where mutation of PatS is reported to lead to a larger number of heterocysts . Table 1 shows that for patS knock-outs in systems with 100–50 cells is lower than 10 cells, and very similar for all systems. In fact, the frequency of observed <10 is around 90.1%, almost twice as much as the wild type with normal patS expression. In addition, standard deviations for patS knock-out are higher than those for the wild type, indicating that knocking-out patS produces an increase in the variability of . Our simulated results for the distribution of L H in patS knock-outs are consistent with the heterocyst spacing distribution reported for a patS deletion strain 48 h after nitrogen step down. In this case the number of vegetative cells between hetercysts ranges from 2 to 14 cells, peaking around 5 cells. Heterocysts spacing distribution becomes broader compared to wild type, peaking towards small intervals. The model however, does not match the exact distribution observed experimentally after prolonged culturing, which becomes broader as the culture progresses  for the same reasons pointed out in the previous case.
Over-expression of patS was also simulated. To our knowledge, this condition has not been explored theoretically with a mathematical model for heterocyst differentiation. In this scenario, simulation results obtained with the present model exhibit a behavior equivalent to that shown in Figure 6-A, i.e. higher expression of patS completely inhibits differentiation into heterocysts.
HetR has been identified as a key factor in heterocyst differentiation. The effect of knocking-out and over-expressing hetR was also simulated, considering an array of 100–50 cells with D = with 5000 randomly chosen initial conditions.
As expected, if hetR is knocked-out no heterocysts are formed . However, hetR over-expression leads to heterocyst differentiation and the formation of heterocyst clusters as shown in Figure 6-D. The pattern obtained is similar to the one observed for patS knock-out. Figure 6-C and 6D respectively, show that the knock-out of patS has a stronger effect than the over-expression of hetR on differentiation into heterocysts, compared to the wild type in Figure 6-B. Results in Figure 6-D are supported by L H histograms shown in Figure 4-C. These results indicate that over-expression of hetR also shifts L H distribution, lowering the average distance between heterocysts , compared to the wild type (see Figure 4-A). This phenomenon has also been observed experimentally [18, 33].
There are no previous reports on the theoretical exploration of the over-expression of hetR using a mathematical model for the gene network involved in heterocyst differentiation. Over-expression of hetR may also lead to the formation of heterocysts under normally repressing conditions when the gene is controlled by an external factor. In the present system, this corresponds to an over-expressed hetR scenario where biologically feasible initial conditions are only restricted to ntcA being off. Simulation results for this case are shown in Figure 6-E: heterocyst formation is observed even in the presence of nitrogen (23.91% less heterocysts than in the absence of nitrogen, Figure 6-D), which is consistent with the experimental observations reported .
Additional factors such as cell growth and death were not considered as they do not affect directly the main activator-inhibitor mechanism. Three types of parameters are considered in this model; parameters related to transport (D), dynamics (matrix A) and a threshold (θ). The system shows to be robust for perturbations in the first two types of parameters, and not surprisingly, very sensitive to changes in θ, which results to be critical as it represents the limit that defines wether a particular element in the network is active or not. It is possible to capture the stable states of the system if the parameters are chosed adequately considering a biological interpretation of the network. The model allows flexibility in the selection of values for the matrix A.
Although the parameters used have not been biologically validated, they retain biological meaning (e.g. >0 for activation, <0 for inhibition and relative magnitudes imply relative effects). In addition, given the characteristics of the implementation, the system can be extended to two and three dimensions.
It is important to note that this model, as other activation-inhibition type models, does not consider any assumptions regarding PatS. It is often said that it is essential that PatS remains inactive until processed in a neighboring vegetative cells for heterocyst-forming cyanobacteria to be able to succesfully activate hetR and therefore achieve their characteristic differentiation pattern. Here, no additional assumptions are considered on PatS, and in fact PatS may act on the cell on which it is produced and still a heterocyst develops, as long as the inhibitory effect is smaller than the counteracting activators effect acting on the same cell.
Additionally, no restrictions are imposed on the range of action of the differentiation activator HetR and the inhibitor PatS. Regardless of these assumptions, the model successfully captures the differentiation pattern observed experimentally for the initial stages of culture. This indicates that it may be possible for the system to achieve the characteristic pattern of heterocyst spacing, without a refractability condition on heterocysts to the action of PatS as discussed in most reviews. Even if PatS inhibited hetR self-activation in the cell where it is produced, pattern formation would be possible if the magnitude of the inhibitory effect on hetR was smaller than hetR self-activation.
One interesting thing is that the pattern formation characteristics naturally emerge from the system without data fitting, since the only parameter that was tuned is the transport factor D. This evolution model was not intended to capture the transient dynamics of the system, but to capture qualitatively the final spacing distribution observed. The simple gene network proposed is able to capture some of the main characteristics of the spacing pattern observed for heterocyst differentiation, with good agreement between simulation and experimental results both for wild type and for perturbations in the gene network.
In this paper the differentiating behavior of cyanobacteria into heterocysts has been studied by means of a mathematical model. A simple gene network which is capable of capturing the complexity of the differentiation process in this cellular system is proposed.
Specific characteristics of the behavior observed experimentally in the presence/absence of nitrogen can be reproduced by the proposed system. The model presented explains how heterocyst-forming cyanobacteria counts to 10, the basic regulatory network elements required and the effect of diffusion on the spacing pattern observed. However, the reason why the system has evolved to this characteristic heterocyst distance remains to be explained. It appears that this has to do with the diffusional characteristics of the inhibitor PatS and the trade-off between resource requirements for heterocysts differentiation vs. their nitrogen fixing capabilities.
In particular, the model reports an average distance between two heterocysts equal to 10 cells, which is one of the most significant characteristics displayed by these systems. In addition, standard deviations of this magnitude were in accordance with the variability observed experimentally [2, 8].
Specific details were also investigated, such as ntcA profiles between two consecutive heterocysts. In this case, simulations results also show a good agreement with experimental observations. Over-expression and knock-out of patS and hetR were also studied. The model is capable of simulating the qualitative behavior of heterocyst-forming cyanobacteria for these scenarios, capturing their distinctive characteristics, for instance: formation of high density clusters of heterocyst or complete inhibition of differentiation in the case of knocked-out and over-expression of patS, respectively. In the case of hetR the opposite behavior is observed: knocking-out hetR completely prevents the formation of heterocysts. On the other hand, over-expression of hetR leads to the formation of heterocyst clusters, but with a lower density and size than observed for patS knock-out. These results validate the proposed network structure.
Extracellular interaction matrix
Intracellular interaction matrix
Global interaction matrix
Cellular state vector
External input vector
Threshold levels vector
Normalized transport factor
This work was supported by the Millennium Scientific Initiative ICM P05-001F. Support deom Fondecyt, Chile grants 1090156 and 1061263 is also acknowledged.
This article has been published as part of BMC Bioinformatics Volume 10 Supplement 6, 2009: European Molecular Biology Network (EMBnet) Conference 2008: 20th Anniversary Celebration. Leading applications and technologies in bioinformatics. The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/10?issue=S6.
This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.