BACKGROUND
the availability of high-throughput techniques for gene expression analysis and identification of promoter-transcription factor  interactions has led to characterization of the intricate gene regulatory networks that govern organism behavior  <cit> .
these networks are composed of a large number of small and topologically distinct subnetworks, including the overrepresented 'network motifs'  <cit> .
in recent years, dynamical systems modeling of regulatory and signaling pathways has provided insight into the equilibrium states and transient dynamics of such subnetworks  <cit> ; for example, detailed cellular and subcellular models demonstrate that interconnected positive and negative feedback loops may give rise to the phenomena of oscillations, excitability, and the existence of multiple stable equilibria   <cit> .
bistability in particular is ubiquitous in biological systems ranging from biochemical networks to ecosystems  <cit> .
in bistable systems, graded inputs  are converted into a discontinuous on/off response  <cit> .
switch-like behavior is also a characteristic of many developmental processes, and it has been suggested that the maintenance of two distinct phenotypic states in the absence of genetic or environmental differences may sometimes be attributed to bistability in an underlying gene network  <cit> .
an intriguing system that exhibits phenotypic bistable behavior, and as such is an excellent candidate for the study of the potential role of bistability in cell fate determination, is the formation of unicellular leaf hairs  in the model plant arabidopsis thaliana.
in arabidopsis, trichomes differentiate from pluripotent epidermal cells by the action of regulatory proteins belonging to the r2r3-myb  and basic helix-loop-helix   classes  <cit> .
these positive regulators directly regulate other tfs  that positively induce trichome initiation, as well as small inhibitory proteins .
a simplified version of this network is shown in fig.
<dig>  an important aspect of trichome differentiation not indicated in fig.
<dig> is the free movement of inhibitory proteins to adjacent epidermal cells  <cit>  where they prevent gl <dig> from interacting with gl <dig>  <cit> , thus creating a domain of surrounding cells that will not become trichomes and resulting in a characteristic spatial pattern.
previous attempts at modeling this cell fate determination system have aimed at explaining how trichome patterns form out of a field of initially equivalent epidermal cells, but have ignored the question of how the primary identity decision is made  <cit> .
such models assume an underlying mechanism of either the 'activator-inhibitor' or 'trapping/depletion' type, both of which include positive regulation of gl <dig> by the gl3-gl <dig> active complex.
consistent with the activator-inhibitor model  <cit> , it has been shown experimentally that the activators do positively control the diffusible inhibitors  <cit> .
however, although recent work has suggested that a positive feedback loop may be involved in cell fate specification in the root epidermis  <cit> , no such direct positive feedback has been found in the leaf.
indeed, the regulators of trichome patterning in the leaf may in fact be involved in a negative feedback loop, as recently demonstrated for gl <dig>  <cit> .
this poses the question of whether the regulatory subnetworks known to be involved in trichome initiation have the capacity for bistability, or whether this is unlikely given the absence of experimental evidence for direct positive autoregulation.
tools for establishing multistability in ode models of regulatory networks
determining if a given gene regulatory network has multiple equilibria first requires the construction of an appropriate mathematical model of the system.
we consider deterministic ode models that take the form  <cit>
where c = t is a vector of concentrations of n species and f = ,⋯,fn)t is a continuous and differentiable matrix function that gives the rate of change of the concentrations over time.
in this context, a biologically meaningful equilibrium solution of eq.
<dig> is a set of  concentrations leading to dc/dt =  <dig>  that is, a vector css that satisfies f =  <dig>  thus, a network with multiple equilibria is one with two or more distinct vectors c for which eq.
<dig> yields dc/dt =  <dig>  a given equilibrium is stable  provided the jacobian matrix,
has n eigenvalues with negative real part when evaluated at the equilibrium concentrations.
since in practice there is poor knowledge of the rate constants and binding affinities that occur in f , we are primarily interested in assessing a network's capacity for multiple equilibria or bistability without specifying these parameters, though all rate constants are assumed to be positive.
a simple parameter-free graphical representation of the system is given in its interaction graph in which each species is represented as a graph node.
a directed positive edge is drawn from species i to j if ∂fj/∂ci >  <dig>  a directed negative edge is drawn if ∂fj/∂ci <  <dig>  no edge is drawn if ∂fi/∂ci =  <dig>
as originally conjectured by thomas  <cit>  and later proven  <cit> , a necessary condition for multiple equilibria  in eq.
<dig> is the presence of a positive circuit in the interaction graph  for at least one set of species concentrations.
thus, the absence of any positive circuits in a network rules out the capacity for bistability, regardless of the parameters chosen for the model.
another necessary condition for multistationarity proposed by kaufman  <cit>  concerns the nuclei of the interaction graph, which are defined as unions of one or more disjoint circuits involving all the vertices of the interaction graph.
multistationarity requires either  the presence of a variable nucleus , or  the presence of two nuclei of opposite signs, where the sign of a nucleus with p positive circuits is p+ <dig>  <cit> .
networks that do not satisfy either of these criteria may also be eliminated in a search for bistability.
when deterministic ode models  are constructed from specified elementary processes under the assumption of mass-action kinetics, their governing equations takes the special form of a chemical reaction network   <cit> ,
where s is a stoichiometric matrix indicating the change in the number of molecules of each species in a particular reaction, and ϕ is a vector of the reaction rates .
while both eq.
<dig> and eq.
<dig> are nonlinear odes, the specific form of crns enables the application of powerful crn-specific network analysis techniques.
one computational method can rule out multiple equilibria in a crn by determining if the associated polynomial function is injective, i.e., it maps distinct arguments to distinct values  <cit> .
a related graphical technique involves analysis of the crn species-reaction  graph  <cit>  .
there is also a suite of computational methods utilizing the so-called deficiency theorems for crns that are implemented in the publicly-available chemical reaction network toolbox  <cit> .
these deficiency theorems give conditions for the existence, multiplicity, and stability of crn steady-states that may often be applied in a parameter-free context; for example, the deficiency zero theorem states that if a crn has certain topological properties, then within each compatibility class , there is exactly one steady-state with strictly positive concentrations, and this steady-state is locally asympotically stable  <cit> .
interested readers may refer to refs.
<cit>  for details on the full range of theorems implemented in the crnt.
this section has briefly introduced five tools that can be used to determine if a given regulatory network has the capacity for multiple equilibria: the thomas conjecture and kauffman's multistationarity conditions based on analysis of the network's interaction graph , a multistationarity theorem based on the structure of the crn sr graph , a computational method that can establish an injective property for the crn's polynomial function , and crn theory as implemented by the chemical reaction network toolbox .
below we use these tools to analyze several small networks consisting of one gene and one gene product , in order to establish both the effectiveness of the various methods and the subnetworks' capacity for multistability.
we then investigate a large set of more complex two-component subnetworks, including the well-studied 'double negative' feedback system.
these small networks are studied without regard for their frequency of occurrence in the larger regulatory machinery of real biological systems, that is, we do not restrict ourselves to overrepresented network motifs.
lastly, we revisit the epidermal differentiation system in arabidopsis as a real-world example of the applicability of these techniques.
RESULTS
bistability in two positive autoregulatory subnetworks: a comparison of methods
consider a small regulatory network consisting of a single tf gene x that is transcribed and translated into protein p, which in turn positively regulates the production of x by binding to one or more independent cis-regulatory elements in its promoter.
such a module, commonly known as a positive autoregulatory motif and shown schematically in fig.
<dig>  is abundant in the transcriptional regulatory networks of eukaryotes .
simple representations such as that shown in fig.
<dig> hide a significant amount of detail; for example, experimental and computational methods commonly used to establish network architecture usually fail to determine whether tfs bind to dna as monomers, dimers, or as part of higher order structures.
however, in the case of positive autoregulation, the actual form in which p binds to the promoter of x—as a monomer  or as a dimer —has significant implications for the possibility of multistable behavior.
in fig.
3a, the motif is specified as a crn that includes basal production of p , degradation of p , reversible binding of p to the promoter of x , and production of p by the tf-gene complex xp .
in fig.
3b there is also basal production of p from x and degradation of p, however two tfs now associate to form a homodimer , and it is this pp dimer that binds to the promoter of x  and produces p via the dimer-gene complex xpp .
from these sets of elementary reactions, which following ref.
<cit>  can also be visualized as wiring diagrams , one can write an unambiguous system of odes that take the form of a crn .
although it has been shown that transcription and elements of post-transcriptional control combine to regulate the protein production rate, and that stochasticity exists at all levels of regulation  <cit> , for simplicity we model protein production as a single deterministic process .
we address the role of translation and mrna degradation in determining the capacity for multiple stable states in the next section.
for both fig.
3a  and fig.
3b , the existence of positive circuits in the interaction graph  means that the ig-t theorem does not preclude multiple equilibria.
neither does the ig-k theorem preclude multiple equilibria; in fig.
3a, there is a variable nucleus that includes the variable circuit x-p-x, and in fig.
3b, there are two nuclei of opposite signs .
and  +   does not change the results of the ig-t and ig-k analysis; see additional file  <dig> ) the srg method is also unable to rule out bistability in these motifs, in fig.
3a due to the splitting of the x + p complex pair  associated with reaction x → x + p by two even-cycles , and in fig.
3b because the cycle containing p and p + p ⇌ pp is not a 1-cycle or an odd-cycle.
in addition, an injective property cannot be established for either model's polynomial function.
unlike the ig-t, ig-k, srg, and inj methods, the crnt determined that the monomer model cannot admit multiple equilibria for any choice of parameters.
conversely, the crnt found that the dimer model does admit multiple equilibria and generated a set of parameters for which the motif is bistable.
the dimer model mass-action odes are shown in fig.
<dig>  along with a bifurcation diagram of equilibrium species concentrations as a function of the dimer association rate constant .
while it has been previously demonstrated through analysis of a positive autoregulatory motif model with explicit mrna concentration that dimerization is necessary for bistability  <cit> , it is significant that modeling transcription and translation as a single operation leads to similar results.
the bistability of the dimer version of the autoregulatory motif, confirmed using crnt, is consistent with the inability of the ig-t, ig-k, srg, and inj methods to rule it out.
however, these methods are also unable to rule out bistability in even the simple case of positive autoregulation by a protein monomer, suggesting that they may be of limited use in the analysis of more complex systems.
analysis of twelve one-component subnetworks
we extended the preceding analysis to ten additional one-component subnetworks constructed from various subsets of the elementary reactions listed in table  <dig> .
table  <dig> shows that the ig-t, ig-k, and srg methods give little information about the capacity for multistability in these twelve simple subnetworks.
indeed, all that can be determined using these methods is the trivial result that multiple equilibria requires the reversible binding of tfs to the promoters of the component genes; constitutive expression alone is not sufficient.
the ig-t and ig-k methods can only rule out mutistability for subnetwork ab .
the srg method only rules out multistability in subnetworks ab and abe since the sr graphs contain only a single cycle and thus there are no split c-pairs .
seven reactions that are combined to generate twelve subnetworks with one gene and one gene product.
analysis of twelve gene regulatory subnetworks that consist of one gene and one gene product.
the letters that make up the subnetwork names refer to the constituent reactions listed in table  <dig>  the entries yes+/- indicate productive/unproductive binding to the gene promoter.
the five analytical tools ig-t, ig-k, srg, inj, and crnt are described in the text.
in the three cases where multiple equilibria were not ruled out , the crnt provided example sets of rate constants leading to bistability.
in contrast to these graphical methods, the inj ruled out multistability in six of twelve cases, and the crnt ruled out multistability in nine of twelve cases.
furthermore, for the three subnetworks where multiple equilibria were not ruled out by the crnt , parameter sets leading to bistability were provided.
in these two methods we begin to see the strength of the full chemical reaction network theory: with the assumption of mass-action kinetics, the underlying structure of crns limits the nonlinearities that appear in the full ode model  <cit> .
the bifurcation diagram for subnetwork abefg is shown in fig.
<dig>  similar results for subnetworks abcdefg and abcefg are shown in additional file  <dig>  fig.
s <dig>
as previously mentioned, the twelve subnetworks presented in table  <dig> were constructed under the assumption that, for the purpose of analyzing subnetwork equilibria, transcription and translation may be combined into a single operation.
we repeated our analysis using one-component subnetworks augmented by elementary reactions that make these processes explicit.
the elementary reactions x → x + p, xp → xp + p, xpp → xpp + p were replaced by reactions representing transcription , translation , and mrna degradation .
application of the crnt to these augmented one-component systems led to no change in the results of table  <dig>  and thus we model protein production as a single one-step process for the remainder of this paper.
survey of two-component regulatory subnetworks
consideration of small networks consisting of two genes  and two gene products  is a natural extension of the analysis presented above.
with a second component included, the number of relevant subnetworks increases to  <dig>  .
the subnetworks may be grouped by their number of dimer/promoter association reactions; every two-component subnetwork has between zero and six such reactions .
for bistability analysis, we applied the inj method in an automated fashion using matlab scripts provided by g. craciun  <cit> .
no attempt was made to apply the ig-t, ig-k, and srg methods that were found to be uninformative with respect to subnetworks with one-component .
approximately  <dig>  subnetworks were randomly selected for study, sampled in a manner that ensured that subnetworks with zero to six dimer-binding reactions were represented in the same proportion as in the complete set of subnetworks.
surprisingly, we found that the percentage of two-component subnetworks for which bistability could be ruled out using the inj method drops to only  <dig> %, with a 95% bootstrap confidence interval , down from 50% in the one-component systems.
as shown in table  <dig>  partitioning the results according to the number of dimer-binding reactions reveals the importance of tf dimerization and feedback to multistability: the group with no dimers binding has the highest percentage of subnetworks that cannot be bistable , and as the number of dimer-binding reactions increases, the percentage of subnetworks definitively without multiple equilibria decreases.
multistability could not be ruled out for any of the sampled two-component subnetworks that contain six dimer-promoter binding reactions.
the percentage of two-component subnetworks for which multiple equilibria can be ruled out, based on inj analysis of ~ <dig>  randomly selected subnetworks .
results are grouped according to the number of dimer-binding reactions in each subnetwork.
brackets indicate bootstrap 95% confidence intervals calculated by sampling with replacement  <dig>  times.
overall, the capacity for multiple equilibria is ruled out by the inj method for  <dig> %  of sampled subnetworks.
in the preceding analysis, we did not distinguish between productive and unproductive tf binding reactions.
we define a productive reaction as one in which expression of the corresponding gene is activated.
similarly, an unproductive reaction is one in which gene expression is repressed.
when the analysis is applied to subsets of systems lacking either productive or unproductive tf binding, the percentage of sampled subnetworks for which inj could note rule out bistability increases as the number of binding reactions increases .
note that the inj method is more successful at ruling out bistability for subnetworks lacking productive tf binding than it is for subnetworks lacking unproductive tf binding, perhaps because bistable networks are enriched in productive tf-promoter interactions.
the crnt was found to be the most effective method for establishing and ruling out bistability in the one-component subnetworks.
however, the current release of the crnt cannot be automated, and is thus difficult to apply to a large set of two-component systems.
for this reason we randomly selected  <dig> two-component subnetworks and applied the inj and crnt methods to compare their effectiveness.
the crnt found  <dig> subnetworks that do not support multiple equilibria , and it generated parameter sets that lead to bistability for  <dig> of the remaining  <dig> subnetworks .
this is quite similar to the fraction of one-component subnetworks that were found by the crnt to be bistable  and well within the confidence interval.
conversely, the inj method was only able to rule out bistability for  <dig> of the  <dig> subnetworks ruled out by the crnt .
although derived from a small sample, taken together these results suggest that the inj method is debilitated by network size to a larger degree than the crnt method.
the five two-component subnetworks that crnt established can support bistability  are compared in fig.
<dig>  it can be seen that all five subnetworks contain binding of the x <dig> gene by p <dig> and dimerization of p <dig> by both p <dig> and p <dig>  the basal production and degradation reactions are also common to all but are not shown for clarity.
a number of other reactions are common to various subnetworks.
it would appear that one common feature between all of these bistable systems is sequestration of the productive tf into inactive complexes.
in subnetwork abcejk for example, p <dig> binds productively to x <dig>  but it also complexes with itself, with p <dig>  and it binds unproductively to x <dig>  p <dig> also binds unproductively to x <dig>  competing with p <dig>  interestingly, only adijkqw includes a tf dimer-promoter binding reaction, shown to be essential for bistability in the one-component systems.
though not comprehensive, this analysis does suggest the existence of a large number of different network architectures capable of supporting bistability.
the canonical reciprocal repression subnetwork
we turn our attention to the 'double-negative' feedback system  first described by monod and jacob  <cit> : a small network consisting of pair of genes  in which each gene's product inhibits transcription of the other gene .
it has been shown both theoretically  <cit>  and experimentally  <cit>  that cooperativity in tf binding is required for this subnetwork to exhibit bistable behavior.
an analysis of crn specifications of this small network is thus an excellent test of the techniques discussed herein.
we constructed two representatives of this subnetwork: one with regulation via tf monomers  and the other with tf dimerization and regulation via tf dimer-promoter interactions .
both models also contained basal production  and degradation  of tfs.
crnt analysis confirmed that the subnetwork containing tf dimerization exhibits bistability for a suitable set of parameters, while the subnetwork with double-negative feedback via tf monomers does not have this capacity.
revisiting the trichome differentiation subnetwork with the crnt
as described above and shown in fig.
<dig>  the trichome differentiation subnetwork consists of regulators belonging to the r2r3-myb  and bhlh  classes that directly regulate gl <dig> and cpc among several other targets  <cit> .
as there is no experimental evidence for the direct positive feedback that has been posited in current models of trichome initiation, the question of whether the system can support bistability without it is one of great interest.
we specify the core trichome differentiation subnetwork as a crn that contains reaction equations for basal expression and degradation of cpc, formation of the gl3-gl <dig> and gl3-cpc complexes, and binding of a single gl3-gl <dig> complex to the promoter of cpc to positively regulate its production.
the constituent equations are given in table  <dig> .
on application of the crnt, we found that this subnetwork specification cannot support bistability.
however, we have shown that details of tf-promoter interactions can be a critical determinant for a regulatory network to exhibit multiple stable equilibria.
we thus constructed three additional crn specifications of the trichome differentiation subnetwork that allow for the possibility of dimerization of the gl3-gl <dig> active complex.
as with the monomer-binding model described above, each of these models contains basal expression and degradation of cpc, and gl3-gl <dig> and gl3-cpc complex formation.
in model  <dig>  gl3-gl <dig> can only bind to the promoter of cpc and activate its transcription as a complex dimer  that forms prior to tf-promoter interaction.
model  <dig> contains sequential binding of two gl3-gl <dig> complexes to the separate sites in the cpc promoter, with cpc protein being produced only when both complexes are bound.
model  <dig> is a hybrid of model  <dig> and model  <dig>  where the gl3-gl <dig> can bind to the cpc promoter as a monomer or a dimer, and both can activate transcription of cpc.
the constituent equations for these models are also given in table  <dig>  in contrast to the monomer-binding case, we found models 2- <dig> are all capable of exhibiting bistability, even in the absence of direct positive autoregulation.
constituent reactions for four specifications of the trichome differentiation network .
CONCLUSIONS
we have used a modeling framework in which small gene regulatory networks are specified as chemical reaction networks , with reactions describing transcription factor  production, degradation, dimerization, and association of both monomer and dimer tfs with promoters in a productive or unproductive fashion.
we used these methods to survey twelve subnetworks consisting of one tf gene and gene product  and found that bistability is only exhibited by those that include positive feedback with cooperative tf binding .
this is not unexpected, as direct positive feedback and non-linearity are often prerequisites for bistability in simple systems  <cit> .
when the analysis tools were applied to a large random sample of two-component subnetworks , we found that the presence of cooperative tf binding makes it more difficult to rule out bistability , although it is ruled out more often in subnetworks lacking direct positive feedback than it is for those lacking direct negative feedback .
interestingly, we found several examples of bistable two-component subnetworks that do not include promotor-tf dimer interactions at all , consistent with prior work demonstrating that sequestration and titration of active tfs into inactive complexes can give rise to non-linearity and bistable behavior  <cit> .
our application of the various equilibria analysis tools suggests that the interaction graph-based techniques  that are generally applicable to deterministic ode systems are rarely informative with respect to small gene regulatory networks.
we also found that the deficiency theorems implemented in the chemical reaction network toolbox  are generally more informative than computational  and graphical  methods that attempt to rule out multiple equilibria by establishing injectivity of the reaction network polynomial function.
we emphasize that the success of the crnt is partly due to the special form of chemical reaction networks that assume kinetics of the mass-action form.
however, this assumption would appear to be a valid one under many working conditions.
that the crnt accurately determined the stability properties of the well-studied positive autoregulatory and reciprocal repression systems , and that one-component systems found by crnt to be bistable were confirmed as such with bifurcation diagrams , is further validation of the crnt's usefulness.
it is important to note that the analytical methods described herein are not the only ones available for probing the stability properties of a dynamical system.
one additional method involves evaluation of the jacobian matrix , because the existence of a saddle-node bifurcation leading to stable and unstable equilibria requires the jacobian to be singular at the critical point css } = 0)  <cit> .
we do not use continuation methods to directly rule out saddle-node bifurcations in this work, as our attempts at such analysis have not yielded results, and to our knowledge it would be difficult to explore large numbers of regulatory networks in this fashion.
furthermore, deterministic odes are only one of several mathematical structures used to model gene regulatory networks; other possibilities include boolean networks, petri nets, markov chains, stochastic odes, and so on.
still, the parameter-free approach to modeling that we have emphasized here is an attractive one for at least two reasons:  there are often multiple assumptions made  when complex network models are constructed as a generic system of odes, and these assumptions may result in inconsistencies  <cit> , and  as we have shown, crn theory facilitates analysis that is usually more fruitful than techniques available for generic odes.
the current state of crn-based modeling and analysis does not include stochastic effects, despite the fact that there is both experimental and theoretical evidence that low concentrations of tfs or binding sites can lead to fluctuations  that may sometimes be an essential aspect of network function  <cit> .
the crn approach thus represents a first step towards generating mixed models of biological systems that incorporate stochastic effects.
future work could extend crn-based techniques to include stochastic dynamics consistent with the elementary processes that compose gene regulatory networks  <cit> .
we further anticipate that new releases of the crnt software package will also be amenable to automation and will allow for crns of unlimited size.
although the modeling formalism used here is a simple one, it is not clear that an increase in complexity would lead to different conclusions; for example, we have already demonstrated that excluding mrna does not influence the one-component subnetwork results.
two other simplifications made find significant experimental support:  that tf dimers are stable to proteolytic degradation , a reasonable assumption given that protein dimerization typically provides protection against proteolysis by enhancing thermal stability or through blocking access to monomer degradation tags  <cit> , and  the exclusion of higher-order oligomeric tf complexes , an assumption appropriate for many eukaryotic tf gene families where dimerization is sufficient for functionality  <cit> .
on the other hand, tf activation can occur by posttranslational modification , a situation analogous to cell signaling pathways in which dual phosphorylation/dephosphorylation cycles confer bistability  <cit> .
while it is possible to extend the modeling formalism illustrated here in these and other directions, the additional combinatorial complexity would make surveys of multi-component networks computationally challenging.
applying the crnt to the real-world example of arabidopsis epidermal cell differentiation, we have shown that gene regulatory subnetworks consistent with established components and interactions may have the capacity for bistability.
this capacity of course depends on the specific details of the network architecture.
in particular, we found that dimerization of the active gl3-gl <dig> complex and cooperative sequential binding of gl3-gl <dig> to the cpc promoter are each independently sufficient to generate a bistable subnetwork, even in the absence of direct positive autoregulation.
both of these alternatives have some experimental support.
sequential binding of gl3-gl <dig> to the cpc promoter is plausible, as multiple tf binding sites in gene promoters are a general characteristic of eukaryotic gene expression  <cit> .
and while dimerization of gl3-gl <dig> has not been experimentally demonstrated, gl <dig> can dimerize  <cit>  either through the conserved bhlh domain, as is the case for other bhlh factors  <cit> , or through the conserved c-terminal act domain  <cit> .
bistability may also be achieved if the gl3-gl <dig> complex is not stable against proteolysis .
if it can be determined that gl3-gl <dig> does not in fact dimerize or bind to two separate binding sites in the cpc promoter, and that it is in fact a stable complex, then the current view of the core trichome differentiation network  must be incomplete.
in that case, there are several intriguing possibilities.
it may be that one or more of the many additional trichome-related genes that have been recently identified  <cit>  are as central to trichome differentiation as gl <dig>  gl <dig>  and cpc.
it is possible that epidermal cell differentiation in arabidopsis is not in fact subserved by a bistable gene regulatory network.
it may also be that some unknown extrinsic mechanism plays a role in arabidopsis epidermal cell fate determination.
a mechanism of this kind may render a bistable network irreversible by causing one equilibrium state to be inaccessible  <cit> , or it may fix a transient activated state before the system can return to a basal equilibrium  <cit> .
there is some evidence for these two latter possibilities, as a differentiated trichome cell is unlikely to revert back to an undifferentiated state.
in particular, a recent study has shown that the induction of gl <dig> function for  <dig> hrs is sufficient to trigger trichome initiation, with development continuing even after gl <dig> function is completely removed  <cit> .
clearly there are a number of possibilities that we have only begun to investigate.
future experiments will be required to distinguish between them and may yet provide insights into this interesting developmental model system.

 the availability of high-throughput techniques for gene expression analysis and identification of promoter-transcription factor interactions has led to characterization of the intricate gene regulatory networks that govern organism behavior. these networks are composed of a large number of small and topologically distinct subnetworks, including the overrepresented'network motifs'. in recent years, dynamical systems modeling of regulatory and signaling pathways has provided insight into the equilibrium states and transient dynamics of such subnetworks; for example, detailed cellular and subcellular models demonstrate that interconnected positive and negative feedback loops may give rise to the phenomena of oscillations, excitability, and the existence of multiple stable equilibria. bistability in particular is ubiquitous in biological systems ranging from biochemical networks to ecosystems. in bistable systems, graded inputs are converted into a discontinuous on/off response. switch-like behavior is also a characteristic of many developmental processes, and it has been suggested that the maintenance of two distinct phenotypic states in the absence of genetic or environmental differences may sometimes be attributed to bistability in an underlying gene network. an intriguing system that exhibits phenotypic bistable behavior, and as such is an excellent candidate for the study of the potential role of bistability in cell fate determination, is the formation of unicellular leaf hairs in the model plant arabidopsis thaliana. in arabidopsis, trichomes differentiate from pluripotent epidermal cells by the action of regulatory proteins belonging to the r2r3-myb and basic helix-loop-helix classes. these positive regulators directly regulate other tfs that positively induce trichome initiation, as well as small inhibitory proteins. a simplified version of this network is shown in fig. an important aspect of trichome differentiation not indicated in fig. is the free movement of inhibitory proteins to adjacent epidermal cells where they prevent gl from interacting with gl, thus creating a domain of surrounding cells that will not become trichomes and resulting in a characteristic spatial pattern. previous attempts at modeling this cell fate determination system have aimed at explaining how trichome patterns form out of a field of initially equivalent epidermal cells, but have ignored the question of how the primary identity decision is made. such models assume an underlying mechanism of either the'activator-inhibitor' or'trapping/depletion' type, both of which include positive regulation of gl by the gl3-gl active complex. consistent with the activator-inhibitor model, it has been shown experimentally that the activators do positively control the diffusible inhibitors. however, although recent work has suggested that a positive feedback loop may be involved in cell fate specification in the root epidermis, no such direct positive feedback has been found in the leaf. indeed, the regulators of trichome patterning in the leaf may in fact be involved in a negative feedback loop, as recently demonstrated for gl. this poses the question of whether the regulatory subnetworks known to be involved in trichome initiation have the capacity for bistability, or whether this is unlikely given the absence of experimental evidence for direct positive autoregulation. tools for establishing multistability in ode models of regulatory networks determining if a given gene regulatory network has multiple equilibria first requires the construction of an appropriate mathematical model of the system. we consider deterministic ode models that take the form where c = t is a vector of concentrations of n species and f =,⋯,fn)t is a continuous and differentiable matrix function that gives the rate of change of the concentrations over time. in this context, a biologically meaningful equilibrium solution of eq. is a set of concentrations leading to dc/dt = that is, a vector css that satisfies f = thus, a network with multiple equilibria is one with two or more distinct vectors c for which eq. yields dc/dt = a given equilibrium is stable provided the jacobian matrix, has n eigenvalues with negative real part when evaluated at the equilibrium concentrations. since in practice there is poor knowledge of the rate constants and binding affinities that occur in f, we are primarily interested in assessing a network's capacity for multiple equilibria or bistability without specifying these parameters, though all rate constants are assumed to be positive. a simple parameter-free graphical representation of the system is given in its interaction graph in which each species is represented as a graph node. a directed positive edge is drawn from species i to j if ∂fj/∂ci > a directed negative edge is drawn if ∂fj/∂ci < no edge is drawn if ∂fi/∂ci = as originally conjectured by thomas and later proven, a necessary condition for multiple equilibria in eq. is the presence of a positive circuit in the interaction graph for at least one set of species concentrations. thus, the absence of any positive circuits in a network rules out the capacity for bistability, regardless of the parameters chosen for the model. another necessary condition for multistationarity proposed by kaufman concerns the nuclei of the interaction graph, which are defined as unions of one or more disjoint circuits involving all the vertices of the interaction graph. multistationarity requires either the presence of a variable nucleus, or the presence of two nuclei of opposite signs, where the sign of a nucleus with p positive circuits is p+. networks that do not satisfy either of these criteria may also be eliminated in a search for bistability. when deterministic ode models are constructed from specified elementary processes under the assumption of mass-action kinetics, their governing equations takes the special form of a chemical reaction network, where s is a stoichiometric matrix indicating the change in the number of molecules of each species in a particular reaction, and ϕ is a vector of the reaction rates. while both eq. and eq. are nonlinear odes, the specific form of crns enables the application of powerful crn-specific network analysis techniques. one computational method can rule out multiple equilibria in a crn by determining if the associated polynomial function is injective, i.e., it maps distinct arguments to distinct values. a related graphical technique involves analysis of the crn species-reaction graph. there is also a suite of computational methods utilizing the so-called deficiency theorems for crns that are implemented in the publicly-available chemical reaction network toolbox. these deficiency theorems give conditions for the existence, multiplicity, and stability of crn steady-states that may often be applied in a parameter-free context; for example, the deficiency zero theorem states that if a crn has certain topological properties, then within each compatibility class, there is exactly one steady-state with strictly positive concentrations, and this steady-state is locally asympotically stable. interested readers may refer to refs. for details on the full range of theorems implemented in the crnt. this section has briefly introduced five tools that can be used to determine if a given regulatory network has the capacity for multiple equilibria: the thomas conjecture and kauffman's multistationarity conditions based on analysis of the network's interaction graph, a multistationarity theorem based on the structure of the crn sr graph, a computational method that can establish an injective property for the crn's polynomial function, and crn theory as implemented by the chemical reaction network toolbox. below we use these tools to analyze several small networks consisting of one gene and one gene product, in order to establish both the effectiveness of the various methods and the subnetworks' capacity for multistability. we then investigate a large set of more complex two-component subnetworks, including the well-studied'double negative' feedback system. these small networks are studied without regard for their frequency of occurrence in the larger regulatory machinery of real biological systems, that is, we do not restrict ourselves to overrepresented network motifs. lastly, we revisit the epidermal differentiation system in arabidopsis as a real-world example of the applicability of these techniques.  bistability in two positive autoregulatory subnetworks: a comparison of methods consider a small regulatory network consisting of a single tf gene x that is transcribed and translated into protein p, which in turn positively regulates the production of x by binding to one or more independent cis-regulatory elements in its promoter. such a module, commonly known as a positive autoregulatory motif and shown schematically in fig. is abundant in the transcriptional regulatory networks of eukaryotes. simple representations such as that shown in fig. hide a significant amount of detail; for example, experimental and computational methods commonly used to establish network architecture usually fail to determine whether tfs bind to dna as monomers, dimers, or as part of higher order structures. however, in the case of positive autoregulation, the actual form in which p binds to the promoter of x—as a monomer or as a dimer —has significant implications for the possibility of multistable behavior. in fig. 3a, the motif is specified as a crn that includes basal production of p, degradation of p, reversible binding of p to the promoter of x, and production of p by the tf-gene complex xp. in fig. 3b there is also basal production of p from x and degradation of p, however two tfs now associate to form a homodimer, and it is this pp dimer that binds to the promoter of x and produces p via the dimer-gene complex xpp. from these sets of elementary reactions, which following ref. can also be visualized as wiring diagrams, one can write an unambiguous system of odes that take the form of a crn. although it has been shown that transcription and elements of post-transcriptional control combine to regulate the protein production rate, and that stochasticity exists at all levels of regulation, for simplicity we model protein production as a single deterministic process. we address the role of translation and mrna degradation in determining the capacity for multiple stable states in the next section. for both fig. 3a and fig. 3b, the existence of positive circuits in the interaction graph means that the ig-t theorem does not preclude multiple equilibria. neither does the ig-k theorem preclude multiple equilibria; in fig. 3a, there is a variable nucleus that includes the variable circuit x-p-x, and in fig. 3b, there are two nuclei of opposite signs. and + does not change the results of the ig-t and ig-k analysis; see additional file ) the srg method is also unable to rule out bistability in these motifs, in fig. 3a due to the splitting of the x + p complex pair associated with reaction x → x + p by two even-cycles, and in fig. 3b because the cycle containing p and p + p ⇌ pp is not a 1-cycle or an odd-cycle. in addition, an injective property cannot be established for either model's polynomial function. unlike the ig-t, ig-k, srg, and inj methods, the crnt determined that the monomer model cannot admit multiple equilibria for any choice of parameters. conversely, the crnt found that the dimer model does admit multiple equilibria and generated a set of parameters for which the motif is bistable. the dimer model mass-action odes are shown in fig. along with a bifurcation diagram of equilibrium species concentrations as a function of the dimer association rate constant. while it has been previously demonstrated through analysis of a positive autoregulatory motif model with explicit mrna concentration that dimerization is necessary for bistability, it is significant that modeling transcription and translation as a single operation leads to similar results. the bistability of the dimer version of the autoregulatory motif, confirmed using crnt, is consistent with the inability of the ig-t, ig-k, srg, and inj methods to rule it out. however, these methods are also unable to rule out bistability in even the simple case of positive autoregulation by a protein monomer, suggesting that they may be of limited use in the analysis of more complex systems. analysis of twelve one-component subnetworks we extended the preceding analysis to ten additional one-component subnetworks constructed from various subsets of the elementary reactions listed in table. table shows that the ig-t, ig-k, and srg methods give little information about the capacity for multistability in these twelve simple subnetworks. indeed, all that can be determined using these methods is the trivial result that multiple equilibria requires the reversible binding of tfs to the promoters of the component genes; constitutive expression alone is not sufficient. the ig-t and ig-k methods can only rule out mutistability for subnetwork ab. the srg method only rules out multistability in subnetworks ab and abe since the sr graphs contain only a single cycle and thus there are no split c-pairs. seven reactions that are combined to generate twelve subnetworks with one gene and one gene product. analysis of twelve gene regulatory subnetworks that consist of one gene and one gene product. the letters that make up the subnetwork names refer to the constituent reactions listed in table the entries yes+/- indicate productive/unproductive binding to the gene promoter. the five analytical tools ig-t, ig-k, srg, inj, and crnt are described in the text. in the three cases where multiple equilibria were not ruled out, the crnt provided example sets of rate constants leading to bistability. in contrast to these graphical methods, the inj ruled out multistability in six of twelve cases, and the crnt ruled out multistability in nine of twelve cases. furthermore, for the three subnetworks where multiple equilibria were not ruled out by the crnt, parameter sets leading to bistability were provided. in these two methods we begin to see the strength of the full chemical reaction network theory: with the assumption of mass-action kinetics, the underlying structure of crns limits the nonlinearities that appear in the full ode model. the bifurcation diagram for subnetwork abefg is shown in fig. similar results for subnetworks abcdefg and abcefg are shown in additional file fig. s as previously mentioned, the twelve subnetworks presented in table were constructed under the assumption that, for the purpose of analyzing subnetwork equilibria, transcription and translation may be combined into a single operation. we repeated our analysis using one-component subnetworks augmented by elementary reactions that make these processes explicit. the elementary reactions x → x + p, xp → xp + p, xpp → xpp + p were replaced by reactions representing transcription, translation, and mrna degradation. application of the crnt to these augmented one-component systems led to no change in the results of table and thus we model protein production as a single one-step process for the remainder of this paper. survey of two-component regulatory subnetworks consideration of small networks consisting of two genes and two gene products is a natural extension of the analysis presented above. with a second component included, the number of relevant subnetworks increases to. the subnetworks may be grouped by their number of dimer/promoter association reactions; every two-component subnetwork has between zero and six such reactions. for bistability analysis, we applied the inj method in an automated fashion using matlab scripts provided by g. craciun. no attempt was made to apply the ig-t, ig-k, and srg methods that were found to be uninformative with respect to subnetworks with one-component. approximately subnetworks were randomly selected for study, sampled in a manner that ensured that subnetworks with zero to six dimer-binding reactions were represented in the same proportion as in the complete set of subnetworks. surprisingly, we found that the percentage of two-component subnetworks for which bistability could be ruled out using the inj method drops to only %, with a 95% bootstrap confidence interval, down from 50% in the one-component systems. as shown in table partitioning the results according to the number of dimer-binding reactions reveals the importance of tf dimerization and feedback to multistability: the group with no dimers binding has the highest percentage of subnetworks that cannot be bistable, and as the number of dimer-binding reactions increases, the percentage of subnetworks definitively without multiple equilibria decreases. multistability could not be ruled out for any of the sampled two-component subnetworks that contain six dimer-promoter binding reactions. the percentage of two-component subnetworks for which multiple equilibria can be ruled out, based on inj analysis of ~ randomly selected subnetworks. results are grouped according to the number of dimer-binding reactions in each subnetwork. brackets indicate bootstrap 95% confidence intervals calculated by sampling with replacement times. overall, the capacity for multiple equilibria is ruled out by the inj method for % of sampled subnetworks. in the preceding analysis, we did not distinguish between productive and unproductive tf binding reactions. we define a productive reaction as one in which expression of the corresponding gene is activated. similarly, an unproductive reaction is one in which gene expression is repressed. when the analysis is applied to subsets of systems lacking either productive or unproductive tf binding, the percentage of sampled subnetworks for which inj could note rule out bistability increases as the number of binding reactions increases. note that the inj method is more successful at ruling out bistability for subnetworks lacking productive tf binding than it is for subnetworks lacking unproductive tf binding, perhaps because bistable networks are enriched in productive tf-promoter interactions. the crnt was found to be the most effective method for establishing and ruling out bistability in the one-component subnetworks. however, the current release of the crnt cannot be automated, and is thus difficult to apply to a large set of two-component systems. for this reason we randomly selected two-component subnetworks and applied the inj and crnt methods to compare their effectiveness. the crnt found subnetworks that do not support multiple equilibria, and it generated parameter sets that lead to bistability for of the remaining subnetworks. this is quite similar to the fraction of one-component subnetworks that were found by the crnt to be bistable and well within the confidence interval. conversely, the inj method was only able to rule out bistability for of the subnetworks ruled out by the crnt. although derived from a small sample, taken together these results suggest that the inj method is debilitated by network size to a larger degree than the crnt method. the five two-component subnetworks that crnt established can support bistability are compared in fig. it can be seen that all five subnetworks contain binding of the x gene by p and dimerization of p by both p and p the basal production and degradation reactions are also common to all but are not shown for clarity. a number of other reactions are common to various subnetworks. it would appear that one common feature between all of these bistable systems is sequestration of the productive tf into inactive complexes. in subnetwork abcejk for example, p binds productively to x but it also complexes with itself, with p and it binds unproductively to x p also binds unproductively to x competing with p interestingly, only adijkqw includes a tf dimer-promoter binding reaction, shown to be essential for bistability in the one-component systems. though not comprehensive, this analysis does suggest the existence of a large number of different network architectures capable of supporting bistability. the canonical reciprocal repression subnetwork we turn our attention to the'double-negative' feedback system first described by monod and jacob: a small network consisting of pair of genes in which each gene's product inhibits transcription of the other gene. it has been shown both theoretically and experimentally that cooperativity in tf binding is required for this subnetwork to exhibit bistable behavior. an analysis of crn specifications of this small network is thus an excellent test of the techniques discussed herein. we constructed two representatives of this subnetwork: one with regulation via tf monomers and the other with tf dimerization and regulation via tf dimer-promoter interactions. both models also contained basal production and degradation of tfs. crnt analysis confirmed that the subnetwork containing tf dimerization exhibits bistability for a suitable set of parameters, while the subnetwork with double-negative feedback via tf monomers does not have this capacity. revisiting the trichome differentiation subnetwork with the crnt as described above and shown in fig. the trichome differentiation subnetwork consists of regulators belonging to the r2r3-myb and bhlh classes that directly regulate gl and cpc among several other targets. as there is no experimental evidence for the direct positive feedback that has been posited in current models of trichome initiation, the question of whether the system can support bistability without it is one of great interest. we specify the core trichome differentiation subnetwork as a crn that contains reaction equations for basal expression and degradation of cpc, formation of the gl3-gl and gl3-cpc complexes, and binding of a single gl3-gl complex to the promoter of cpc to positively regulate its production. the constituent equations are given in table. on application of the crnt, we found that this subnetwork specification cannot support bistability. however, we have shown that details of tf-promoter interactions can be a critical determinant for a regulatory network to exhibit multiple stable equilibria. we thus constructed three additional crn specifications of the trichome differentiation subnetwork that allow for the possibility of dimerization of the gl3-gl active complex. as with the monomer-binding model described above, each of these models contains basal expression and degradation of cpc, and gl3-gl and gl3-cpc complex formation. in model gl3-gl can only bind to the promoter of cpc and activate its transcription as a complex dimer that forms prior to tf-promoter interaction. model contains sequential binding of two gl3-gl complexes to the separate sites in the cpc promoter, with cpc protein being produced only when both complexes are bound. model is a hybrid of model and model where the gl3-gl can bind to the cpc promoter as a monomer or a dimer, and both can activate transcription of cpc. the constituent equations for these models are also given in table in contrast to the monomer-binding case, we found models 2- are all capable of exhibiting bistability, even in the absence of direct positive autoregulation. constituent reactions for four specifications of the trichome differentiation network.  we have used a modeling framework in which small gene regulatory networks are specified as chemical reaction networks, with reactions describing transcription factor production, degradation, dimerization, and association of both monomer and dimer tfs with promoters in a productive or unproductive fashion. we used these methods to survey twelve subnetworks consisting of one tf gene and gene product and found that bistability is only exhibited by those that include positive feedback with cooperative tf binding. this is not unexpected, as direct positive feedback and non-linearity are often prerequisites for bistability in simple systems. when the analysis tools were applied to a large random sample of two-component subnetworks, we found that the presence of cooperative tf binding makes it more difficult to rule out bistability, although it is ruled out more often in subnetworks lacking direct positive feedback than it is for those lacking direct negative feedback. interestingly, we found several examples of bistable two-component subnetworks that do not include promotor-tf dimer interactions at all, consistent with prior work demonstrating that sequestration and titration of active tfs into inactive complexes can give rise to non-linearity and bistable behavior. our application of the various equilibria analysis tools suggests that the interaction graph-based techniques that are generally applicable to deterministic ode systems are rarely informative with respect to small gene regulatory networks. we also found that the deficiency theorems implemented in the chemical reaction network toolbox are generally more informative than computational and graphical methods that attempt to rule out multiple equilibria by establishing injectivity of the reaction network polynomial function. we emphasize that the success of the crnt is partly due to the special form of chemical reaction networks that assume kinetics of the mass-action form. however, this assumption would appear to be a valid one under many working conditions. that the crnt accurately determined the stability properties of the well-studied positive autoregulatory and reciprocal repression systems, and that one-component systems found by crnt to be bistable were confirmed as such with bifurcation diagrams, is further validation of the crnt's usefulness. it is important to note that the analytical methods described herein are not the only ones available for probing the stability properties of a dynamical system. one additional method involves evaluation of the jacobian matrix, because the existence of a saddle-node bifurcation leading to stable and unstable equilibria requires the jacobian to be singular at the critical point css } = 0). we do not use continuation methods to directly rule out saddle-node bifurcations in this work, as our attempts at such analysis have not yielded results, and to our knowledge it would be difficult to explore large numbers of regulatory networks in this fashion. furthermore, deterministic odes are only one of several mathematical structures used to model gene regulatory networks; other possibilities include boolean networks, petri nets, markov chains, stochastic odes, and so on. still, the parameter-free approach to modeling that we have emphasized here is an attractive one for at least two reasons: there are often multiple assumptions made when complex network models are constructed as a generic system of odes, and these assumptions may result in inconsistencies, and as we have shown, crn theory facilitates analysis that is usually more fruitful than techniques available for generic odes. the current state of crn-based modeling and analysis does not include stochastic effects, despite the fact that there is both experimental and theoretical evidence that low concentrations of tfs or binding sites can lead to fluctuations that may sometimes be an essential aspect of network function. the crn approach thus represents a first step towards generating mixed models of biological systems that incorporate stochastic effects. future work could extend crn-based techniques to include stochastic dynamics consistent with the elementary processes that compose gene regulatory networks. we further anticipate that new releases of the crnt software package will also be amenable to automation and will allow for crns of unlimited size. although the modeling formalism used here is a simple one, it is not clear that an increase in complexity would lead to different conclusions; for example, we have already demonstrated that excluding mrna does not influence the one-component subnetwork results. two other simplifications made find significant experimental support: that tf dimers are stable to proteolytic degradation, a reasonable assumption given that protein dimerization typically provides protection against proteolysis by enhancing thermal stability or through blocking access to monomer degradation tags, and the exclusion of higher-order oligomeric tf complexes, an assumption appropriate for many eukaryotic tf gene families where dimerization is sufficient for functionality. on the other hand, tf activation can occur by posttranslational modification, a situation analogous to cell signaling pathways in which dual phosphorylation/dephosphorylation cycles confer bistability. while it is possible to extend the modeling formalism illustrated here in these and other directions, the additional combinatorial complexity would make surveys of multi-component networks computationally challenging. applying the crnt to the real-world example of arabidopsis epidermal cell differentiation, we have shown that gene regulatory subnetworks consistent with established components and interactions may have the capacity for bistability. this capacity of course depends on the specific details of the network architecture. in particular, we found that dimerization of the active gl3-gl complex and cooperative sequential binding of gl3-gl to the cpc promoter are each independently sufficient to generate a bistable subnetwork, even in the absence of direct positive autoregulation. both of these alternatives have some experimental support. sequential binding of gl3-gl to the cpc promoter is plausible, as multiple tf binding sites in gene promoters are a general characteristic of eukaryotic gene expression. and while dimerization of gl3-gl has not been experimentally demonstrated, gl can dimerize either through the conserved bhlh domain, as is the case for other bhlh factors, or through the conserved c-terminal act domain. bistability may also be achieved if the gl3-gl complex is not stable against proteolysis. if it can be determined that gl3-gl does not in fact dimerize or bind to two separate binding sites in the cpc promoter, and that it is in fact a stable complex, then the current view of the core trichome differentiation network must be incomplete. in that case, there are several intriguing possibilities. it may be that one or more of the many additional trichome-related genes that have been recently identified are as central to trichome differentiation as gl gl and cpc. it is possible that epidermal cell differentiation in arabidopsis is not in fact subserved by a bistable gene regulatory network. it may also be that some unknown extrinsic mechanism plays a role in arabidopsis epidermal cell fate determination. a mechanism of this kind may render a bistable network irreversible by causing one equilibrium state to be inaccessible, or it may fix a transient activated state before the system can return to a basal equilibrium. there is some evidence for these two latter possibilities, as a differentiated trichome cell is unlikely to revert back to an undifferentiated state. in particular, a recent study has shown that the induction of gl function for hrs is sufficient to trigger trichome initiation, with development continuing even after gl function is completely removed. clearly there are a number of possibilities that we have only begun to investigate. future experiments will be required to distinguish between them and may yet provide insights into this interesting developmental model system. 
in the specific case of epidermal differentiation in arabidopsis, dimerization of the gl3-gl complex and cooperative sequential binding of gl3-gl to the cpc promoter are each independently sufficient for bistability. we find that positive feedback and cooperativity mediated by transcription factor dimerization is a requirement for one-component subnetwork bistability. in particular, there has been considerable interest in using mathematical tools to understand how multistable regulatory networks may contribute to developmental processes such as cell fate determination. indeed, such a network may subserve the formation of unicellular leaf hairs in the model plant arabidopsis thaliana. these methods are compared and applied to a full set of one-component subnetworks, as well as a large random sample from similarly constructed two-component subnetworks. for subnetworks with two components, the presence of these processes increases the probability that a randomly sampled subnetwork will exhibit multiple equilibria, although we find several examples of bistable two-component subnetworks that do not involve cooperative tf-promoter binding. in order to investigate the capacity of small gene regulatory networks to generate multiple equilibria, we present a chemical reaction network -based modeling formalism and describe a number of methods for crn analysis in a parameter-free context. computational methods utilizing crn-specific theorems to rule out bistability in small gene regulatory networks are far superior to techniques generally applicable to deterministic ode systems. using these methods to conduct an unbiased survey of parameter-free deterministic models of small networks, and the arabidopsis epidermal cell differentiation subnetwork in particular, we illustrate how future experimental research may be guided by network structure analysis. recent years have seen a dramatic increase in the use of mathematical modeling to gain insight into gene regulatory network behavior across many different organisms.
 recent years have seen a dramatic increase in the use of mathematical modeling to gain insight into gene regulatory network behavior across many different organisms. in particular, there has been considerable interest in using mathematical tools to understand how multistable regulatory networks may contribute to developmental processes such as cell fate determination. indeed, such a network may subserve the formation of unicellular leaf hairs in the model plant arabidopsis thaliana.  in order to investigate the capacity of small gene regulatory networks to generate multiple equilibria, we present a chemical reaction network -based modeling formalism and describe a number of methods for crn analysis in a parameter-free context. these methods are compared and applied to a full set of one-component subnetworks, as well as a large random sample from similarly constructed two-component subnetworks. we find that positive feedback and cooperativity mediated by transcription factor dimerization is a requirement for one-component subnetwork bistability. for subnetworks with two components, the presence of these processes increases the probability that a randomly sampled subnetwork will exhibit multiple equilibria, although we find several examples of bistable two-component subnetworks that do not involve cooperative tf-promoter binding. in the specific case of epidermal differentiation in arabidopsis, dimerization of the gl3-gl complex and cooperative sequential binding of gl3-gl to the cpc promoter are each independently sufficient for bistability.  computational methods utilizing crn-specific theorems to rule out bistability in small gene regulatory networks are far superior to techniques generally applicable to deterministic ode systems. using these methods to conduct an unbiased survey of parameter-free deterministic models of small networks, and the arabidopsis epidermal cell differentiation subnetwork in particular, we illustrate how future experimental research may be guided by network structure analysis. 