article: Heterocyst differentiation in cyanobacteria filaments is one of the simplest examples of cellular differentiation and pattern formation in multicellular organisms. 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 cyanobacterium 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.
  The formation of multicellular organisms from the assembly of single-celled ones constitutes one of the most striking and complex problems tackled by biology. 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). A landmark process of (prokaryotic) cellular differentiation and cooperative pattern formation is the heterocyst differentiation in cyanobacteria filaments [3,4]. 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, nitrogenase, 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 (i) 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-fixing cells distributed along the filament and acting as nitrogen suppliers. In the presence of combined nitrogen (such as nitrate, nitrite, ammonium or urea), most cyanobacteria (Anabaena PCC strain 7120 being the most representative example) form long filaments of photosynthetic vegetative cells. 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]. Let us briefly review the previous studies on the mathematical modeling of heterocyst pattern formation. 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. Our model connects the diffusion of combined nitrogen along the filament with the dynamical properties of the underlying genetic circuit of each single cyanobacterium, capturing both the development of heterocyst patterns and their maintenance. 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. The work is structured as follows. 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.
  Heterocyst development begins with sensing combined-nitrogen (cN) limitation and ends with nitrogen fixation in mature heterocysts. 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. The process is initiated with the accumulation of 2-oxoglutarate (2-OG) as a consequence of cN deprivation [9,16]. 2-OG interacts with ammonium through the GS/GOGAT cycle [17–19] (see Fig. 2). Under cN starvation, the GS/GOGAT cycle breaks down, leading to the accumulation of 2-OG inside the cell [9]. In its turn, 2-OG stimulates the DNA-binding activity of NtcA, an important transcription factor for heterocyst development [18,20,21]. Furthermore, the transcription of the genes targeted by NtcA does not start in the absence of 2-OG [22,23]. NtcA presents autoregulation [22,24,25] and indirectly activates the key gene that controls cell differentiation and pattern formation: hetR [26–28]. To bind DNA, NtcA needs to homodimerize [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]. The next step in heterocyst development is the activation of hetR. 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 homodimerization [34,35]. Multiple transcription factors related to heterocyst formation are up-regulated by HetR, including hetR itself [35], ntcA [36] and patS [35]. The up-regulatory loop of NtcA and HetR is essential for heterocyst differentiation [36–39]. 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 upregulation by HetR [8]. 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 [40]. 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]. The last stages of heterocyst development cause the physiological changes of the cell aimed at creating an anaerobic environment that sustains nitrogen fixation. 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, nif genes, are expressed. These genes encode, among others, the enzyme nitrogenase, which ultimately performs nitrogen fixation. The fixed nitrogen of the new heterocysts acts as an inhibitor of the differentiation together with PatS and HetN [45]. 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 heterocysts produce fixed 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. In this section we translate the genetic circuit previously described into a set of differential equations, for which we follow the derivation in [46–48]. Details are left to supplementary information (S1 text). To simplify notation, constants related to NtcA, HetR, PatS and cN are denoted with the letters a, r, s, and n respectively. We start by looking at the transcription of ntcA, which is regulated by HetR and NtcA. 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: v a = L a + v a a κ a a [ 2-OG ] [ NtcA ] 2 + v a r κ a r [ HetR ] 2 + v a a r κ a a κ a r [ 2-OG ] [ NtcA ] 2 [ HetR ] 2 (1 + κ a a [ 2-OG ] [ NtcA ] 2) (1 + κ a r [ HetR ] 2), (1) where va measures the production rate of NtcA in units of concentration per time, v a a, v r a and v a r a are the rates when only NtcA, only HetR or both are bound to DNA respectively, and κ * a are the inverse of the effective dissociation constants of the compounds that bind DNA. 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. Similarly we can obtain the transcription velocity for HetR. 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 auto-regulatory 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: v r = L r + v r a κ r a [ 2-OG ] [ NtcA ] 2 (1 + κ r s [ PatS ]) + v r r κ r r [ HetR ] 2 + v r a r κ r a κ r r [ 2-OG ] [ NtcA ] 2 [ HetR ] 2 1 + κ r a [ 2-OG ] [ NtcA ] 2 1 + κ r r [ HetR ] 2 + κ r s [ PatS ] (2) HetR regulates most processes of the genetic circuit. It governs, among others, the transcription of ntcA, patS, hep, hgl and nif genes that lead to most of structural changes of the cell and to nitrogen fixation. The inhibitor PatS is regulated by HetR, and we assume no other influence. 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: v s = L s + v s r κ s r [ HetR ] 2 1 + κ s r [ HetR ] 2 (3) Finally, we have to relate nitrogenase concentration [Ni] to that of combined Nitrogen [cN], both regulated by HetR and the levels of 2-OG [2-OG]. Let us begin by examining nitrogenase concentration, which is directly controlled by nif genes. Although this is not a direct process, we can assume, as we did for the NtcA-regulation of hetR, that nif genes are functionally governed by [HetR] following a typical Hill function. The nitrogenase production rate is given by: d [Ni] d t = L N i + v N i r κ N i r [HetR] 2 1 + κ N i r [HetR] 2 − δ N i [Ni]. (4) where δNi represents the degradation rate of nitrogenase. We can effectively account for the lag introduced by intermediate processes not taken into account explicitly in the model by increasing the value of δNi 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: d [cN] d t = L n ′ + v n ′ [Ni] − δ n ′ [cN], (5) where L n ′ represents the flux of cN from the exterior of the cell. Assuming that the levels of cN relax rapidly we solve Eq. (5) for the steady state. Substituting in Eq. (4) we find: d [cN] d t = L n + v n r κ n r [HetR] 2 1 + κ n r [HetR] 2 − δ n [cN], (6) where Ln=1δn′ (vn′LNi+δNiLn′), vnr=vn′δn′vNir, δn=δNi, κnr=κNir. (7) To get a closed system of equations, we shall investigate the relation between cN and 2-OG. 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: [ glutamate ] = κ ← [ glutamine ] [ 2-OG ], [ glutamine ] = κ → [ glutamate ] [ cN ], (8) which lead to the relation: [2-OG] = 1 κ ← κ → [cN]. (9) However, this expression does not behave properly for small concentrations of cN, which are expected under cN deprivation: 2-OG levels would increase without limit. 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 [cN] in Eq. (9) [2-OG] = 1 κ 2-OG + κ ← κ → [cN], (10) which reaches the maximum value [2-OG]max = 1/κ2-OG at [cN] = 0. Finally we introduce the differential equations governing cyanobacterial reaction to nitrogen deprivation. 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, δ*, we find: dqadτ=la+βaaγaaqa2+βarγarqr2 (1+qn) +βaarγaaqa2γarqr2 (1+qn+γaaqa2) (1+γarqr2) −daqa, dqrdτ=lr+βraqa2 (1+qs) +βrrqr2 (1+qn) +βrarqa2qr2 (1+qn+qa2) (1+qs+qr2) −qr, dqsdτ=ls+βsrγsrqr21+γsrqr2−dsqs, dqndτ=ln+βnrγnrqr21+γnrqr2−dnqn, (11) where we have introduced the dimensionless variables: qa=κraκ2-OG︸ϕa[NtcA], qr=κrr︸ϕr[HetR], qs=κrs︸ϕs[PatS], qn=κ←κ→κ2-OG︸ϕn[cN], τ=δrt, (12) and the constants l * = L * ϕ * δ r, β * • = v * • ϕ * δ r, γ * • = κ * • κ r •, d * = δ * δ r. (13) Let us finally stress that this is a deterministic model for a single cyanobacterium. 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. 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 NctA in their instantaneous equilibrium. 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 qn f a (q s, q n) = d q a d τ = 0, f r (q s, q n) = d q r d τ = 0. (14) The numerical solution to this problem is sketched in Fig. 4. 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. 4I) 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. 4II). In the bistable region the system behaves as a switch that can be either OFF in a vegetative state (lower branch, with a small production of HetR and NtcA) or ON in a heterocyst state (upper branch, with a high production of HetR and NtcA). 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. With (14) solved, we can apply the solution to calculate the effective field sensed by the (qs, qn) pair. 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. In the presence of cN (Fig. 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). Diffusion protects cells in the neighborhood of the newly formed heterocyst to initiate the differentiation: as heterocysts are producers of cN and PatS, the vegetative fixed point of the cells in its neighborhood will move towards an A-like state, thus becoming more stable to perturbations. The heterocyst fixed point also becomes more stable due to diffusion, since its production of inhibitors is distributed among other cells (see Fig. 5D). In the previous section, we introduced a single cell model for the cyanobacteria reaction to nitrogen-limiting conditions. 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. In this section, we extend the previous results and consider a chain of vegetative cells facing nitrogen deprivation. 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: d C i d t = D C C i + 1 + C i − 1 − 2 C i. (15) where DC is called the diffusion 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: dqi, adτ=la+βaaγaaqa2+βarγarqr2 (1+qn) +βaarγaaqa2γarqr2 (1+qn+γaaqa2) (1+γarqr2) −daqi, a+Gi, a (t), dqi, rdτ=lr+βraqa2 (1+qs) +βrrqr2 (1+qn) +βrarqa2qr2 (1+qn+qa2) (1+qs+qr2) −qi, r+Gi, r (t), dqi, sdτ=ls+βsrγsrqr21+Ksrqr2−dsqi, s+Ds (qi+1, s+qi−1, s−2qi, s) +Gi, s (t), dqi, ndτ=ln+βnrγnrqr21+γnrqr2−dnqi, n+Dn (qi+1, n+qi−1, n−2qi, n) +Gi, n (t), (16) which constitutes the model for a cyanobacteria filament. To account for environment variability we add white noise, Gi, * (t), of the same amplitude, ⟨Gi, * (t) Gi, * (t′) ⟩ = ξδ (t − t′), 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 S2 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 (see Methods) [54]. 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 ξ = 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 Ds = 0. 1 and Dn = 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. In Fig. 6 we show the dynamics that the 4 variables exhibit when the filament is deprived of cN. 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. Fig. 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 (i. e. those vegetative cells that are not supplied of sufficient cN and PatS). Finally, a semiregular 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. Finally, in Fig. 7B we show the time-evolution of the histogram for the distance between two consecutive heterocysts. It should be stressed that although initially some close heterocysts appear, they are eliminated by the non-linear 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 Γ-distribution, implying Poisson-distributed wait times for the underlying noise driven process.
  The study of cell differentiation and its underlying mechanisms constitutes one of most intriguing problems in biology. This phenomenon is the basis for multicellular organisms 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 attractors of the non-linear dynamics of the regulatory equations. However, the study of many coupled cells is needed as cyanobacteria do not differentiate when isolated. The model is rounded off by coupling a number of cells in a one-dimensional array so that combined nitrogen and PatS can diffuse along the cellular chain. 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]. Although the model presented here integrates both the internal cell dynamics and the coupling between cells via diffusion, there exist other ingredients that can be also incorporated. 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 Γ-function shape of Fig. 7B and thus approaching better to the experimental distribution. Other improvements to the approach presented here will come from the availability of more experimental data. 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].
  To reproduce the dynamics of Eq. (16) we make use of the integration scheme proposed in [54]. 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 x ˙ = f (x) + G (t), (17) where G (t) is a Gaussian white noise with ⟨ G (t) ⟩ = 0, and, ⟨ G (t) G (t ′) ⟩ = ξ δ (t ′ − t), 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 ξ = 0. In this work we have employed a 3O4S2G algorithm, which is correct up to 3th order, is developed in 4 stages and uses 2 independent Gaussian random variables.
summary: Cyanobacteria filaments are paradigmatic examples of prokaryotic cellular differentiation and cooperative pattern formation. 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.
title: An Integrative Approach for Modeling and Simulation of Heterocyst Pattern Formation in Cyanobacteria Filaments