- Research article
- Open Access

# Analysis of a dynamic model of guard cell signaling reveals the stability of signal propagation

- Xiao Gan
^{1}and - Réka Albert
^{1}Email author

**10**:78

https://doi.org/10.1186/s12918-016-0327-7

© The Author(s). 2016

**Received:**10 June 2016**Accepted:**11 August 2016**Published:**19 August 2016

## Abstract

### Background

Analyzing the long-term behaviors (attractors) of dynamic models of biological systems can provide valuable insight into biological phenotypes and their stability. In this paper we identify the allowed long-term behaviors of a multi-level, 70-node dynamic model of the stomatal opening process in plants.

### Results

We start by reducing the model’s huge state space. We first reduce unregulated nodes and simple mediator nodes, then simplify the regulatory functions of selected nodes while keeping the model consistent with experimental observations. We perform attractor analysis on the resulting 32-node reduced model by two methods: 1. converting it into a Boolean model, then applying two attractor-finding algorithms; 2. theoretical analysis of the regulatory functions. We further demonstrate the robustness of signal propagation by showing that a large percentage of single-node knockouts does not affect the stomatal opening level.

### Conclusions

Combining both methods with analysis of perturbation scenarios, we conclude that all nodes except two in the reduced model have a single attractor; and only two nodes can admit oscillations. The multistability or oscillations of these four nodes do not affect the stomatal opening level in any situation. This conclusion applies to the original model as well in all the biologically meaningful cases. In addition, the stomatal opening level is resilient against single-node knockouts. Thus, we conclude that the complex structure of this signal transduction network provides multiple information propagation pathways while not allowing extensive multistability or oscillations, resulting in robust signal propagation. Our innovative combination of methods offers a promising way to analyze multi-level models.

## Keywords

- Network model
- Discrete dynamic model
- Biological network
- Signal transduction
- Plant signaling
- Attractor
- Stomatal opening
- Network reduction
- Boolean conversion
- Stable motif

## Background

Modeling offers a comprehensive way to understand biological processes by integrating the components involved in them and the interactions between components. Models can recapitulate and explain the emergent outcome(s) of the process [1, 2]. Representing cellular processes that involve many proteins and small molecules by a signal transduction network can reveal indirect relationships between components and provide new insight [3–5]. Such network usually consists of nodes representing biological entities, and edges representing interactions. Once a network has been constructed, dynamic modeling, where each node in the network is associated with a variable representing its abundance or activity, can further describe the behavior of the network. Dynamic models can have continuous variables whose change is described by differential equations [6], discrete variables described by discrete (logical) regulatory functions [7, 8], or a combination of continuous and discrete variables [9]. The major advantage of discrete dynamic and continuous-discrete hybrid models is that they use many fewer parameters than continuous models and thus need less parameter estimation [10–12]. Modeling allows one to analyze the biological system represented by the network *in silico*, when performing the relevant experiment is infeasible. It also helps identify general principles of biological systems [13, 14].

The biological process of stomatal opening in plants is a good example of a complex system wherein modeling leads to significant gain in understanding [15, 16]. Stomata are pores on leaf surfaces that allow the plant to exchange carbon dioxide (CO_{2}) and oxygen with the atmosphere. Stomata are formed by two guard cells that can change shape: swelling of guard cells leads to stomatal opening; their shrinking leads to stomatal closure. The shape of each guard cell is directly controlled by water flow through the membrane, which is in turn controlled by ion flow. Different signals can affect the guard cell, changing its ion concentration in direct and indirect ways, resulting in stomatal opening or closure [17–19]. These signals include light of different wavelengths, CO_{2} concentration in the air, and plant hormones like abscisic acid (ABA). The regulation of stomatal opening is essential to plants, as it controls vital activities like the uptake of CO_{2} for photosynthesis, and the unavoidable water loss through evaporation [20]. Through extensive experimentation over several decades, more than 70 proteins and small molecules have been identified to participate in this process.

_{2}, and ABA. The more than 150 edges are directed and signed, with arrowheads indicating activation and terminal black circles indicating inhibition.

^{31}states. The logical regulatory functions, describing each node’s future state based on the states of the node’s regulators, use a combination of Boolean logic operators (

*And, Or, Not*), algebraic operations, and input-output tables. For example, the regulatory function of PRSL1 is:

Here for simplicity the node states are denoted by the node names; the asterisk in “PRSL1*” indicates that this will be the next state of the PRSL1. The “Or” Boolean operator expresses that either of the blue light receptors, i.e. the phot1 complex or phot2, can independently activate PRSL1.

The Sun et al. model starts from an initial condition representative of closed stomata. Then a combination of the four input signals is applied. Red light, blue light, and ABA are represented as binary variables, and external CO_{2} is represented with three states: 0 (CO_{2} free air), 1 (ambient CO_{2}) and 2 (high CO_{2}). The system’s response is simulated through repetitive re-evaluation of each node’s state until a stable value of stomatal opening is observed. The model successfully captures stomatal opening in response to combinations of the signals. It also successfully reproduces stomatal opening under most of the experimentally studied perturbation scenarios (i.e. genetic knockouts or external supply of components). In total, the model is consistent with 63 out of 66 experimental observations collected by Sun et al. [15]. The model predicts the outcome of a large number of scenarios that have not been explored experimentally so far. It also revealed a gap of knowledge regarding the cross-talk of red light and ABA signaling, and filled it with a newly predicted interaction.

Although the Sun et al. model recapitulates existing knowledge and offers new predictions, the model’s full dynamic repertoire could not be characterized due to its large state space. Instead, Sun et al. focused on tracking the output node, stomatal opening, and a few selected internal nodes, in time. In this paper we apply multiple methods to analyze the model and aim to fully map all its potential long-term behaviors, or in other words, attractors.

## Methods

### Attractors of a dynamical system

An attractor is a set of states from which only states in the same set can be reached. Attractors that consist of a single state are called stable steady states or fixed points; attractors that contain multiple states are called complex attractors or oscillations [10]. In biological networks, attractors often have significant biological meaning. In a cell signaling network, attractors correspond to cell types, cell fates or behaviors [21]. For example, one attractor can represent a healthy differentiated cell, while another attractor can represent an abnormally motile cancer cell [22].

### Update scheme of a discrete time model

In the Sun et al. model, as in most discrete dynamic models, time is an implicit variable. As there is very little information about the kinetics of the nodes in the stomatal opening network, the model incorporates an element of stochasticity in timing. The timing does not affect a system’s fixed point attractors, but it can change the complex attractors and the possibility of reaching a given attractor from a given initial state [10]. In the Sun et al. model, a random–order asynchronous update is used. Specifically, at each time step, a random order of nodes (excluding the four input nodes and the output node stomatal opening) is generated, and each node’s state is reevaluated in this order; stomatal opening is always updated last. In the next time step a different order is generated randomly. In this paper, we use a different type of stochastic update, called general asynchronous update, wherein a randomly selected node is updated at each time step. This is required by the network reduction method we use. Although this theoretically could cause a difference in complex attractors, we will show that in this specific model the two update methods yield the same attractors.

### Network reduction

To reduce the Sun et al. model’s state space, we apply a network reduction method developed by Saadatpour et al. [23] that is proven to preserve the attractors of a Boolean model. Two types of nodes can be reduced (eliminated or merged): source nodes with no incoming edges, and simple mediator nodes that have one incoming and/or one outgoing edge. In the reduction, the source node’s state is directly plugged into the regulatory function of all of its direct successor nodes; then the source node is eliminated. For a simple mediator node with one predecessor (regulator) and one successor (target), its regulator is connected to its target and the mediator node is merged into the regulator. If there is one regulator and several targets of the mediator node, but no direct edges between the regulator and any of the targets, the mediator node is merged into the regulator. Conversely, if there are several regulators and one target of the mediator node, but no direct edges among any of the regulators and the target, the mediator node is merged into its target. Although this method is not proven in the multi-level case, we conjecture that attractors are also conserved for a multi-level model, and will show from the results that in the Sun et al. model this reduction method preserved all attractors.

### Elimination of redundant edges

The italicized words “And”, “Or” and “Not” are Boolean logic operators; the non-italicized words represent node names. In this regulatory function every node is Boolean (binary). The first clause “**NADPH**
And
**AtrbohD/F**” and the second “**NADPH**
And
**AtrbohD/F**
And
**CDPK**” are connected with an “Or” rule, with the result that the node “CDPK” does not have any influence on the outcome. Therefore, we can prune the edge from CDPK to ROS without changing the model’s dynamics. We similarly prune three additional redundant edges.

### Converting a multi-level model to Boolean

There are several possibilities to convert a multi-level model to Boolean [24]. The standard method used in the case of logical models of regulatory networks is the Van Ham mapping [25, 26]. It preserves the dynamics of the original model if the variables in the original model can be represented by integers and if the original model only allows state transitions in which one node changes its state by one level [26]. The Sun et al. model does not satisfy these criteria. However there still is a conclusion that we can use: All types of conversions maintain the fixed points and the reachability of states (i.e. if there is a sequence of state transitions from state A to state B before conversion, there must be a sequence of state transitions from the corresponding state A’ to state B’ after the conversion) [26]. So the worst distortion of attractors due to the conversion is the merging of two complex attractors into one. In this light we choose to use an economic mapping of each multi-level node into as many Boolean nodes as necessary for the binary representation of the corresponding integer. We will show that in this specific model, the conversion did not change the attractors.

### Abbreviations

Full names of the network components denoted by abbreviated node names in Fig. 1

Abbreviation | Full name | Abbreviation | Full name |
---|---|---|---|

14-3-3 protein | 14-3-3 protein that binds to the H | 14-3-3 protein | 14-3-3 protein that binds to phototropin 1 |

ABA | abscisic acid | ABI1 | 2C-type protein phosphatase |

acid. of apoplast | the acidification of the apoplast | AnionCh | anion efflux channels at the plasma membrane |

AtABCB14 | ABC transporter gene AtABCB14 | Atnoa1 | protein nitric oxide-associated 1 |

AtrbohD/F | NADPH oxidase D/F | AtSTP1 | H-monosaccharide symporter gene AtSTP1 |

Ca | Ca | CaIC | inward Ca |

CaR | Ca | carbon fixation | light-independent reactions of photosynthesis |

CDPK | Ca | CHL1 | dual-affinity nitrate transporter gene AtNRT1.1 |

C | intercellular CO | FFA | free fatty acids |

H | the phosphorylated H | H | 14-3-3 protein bound H |

KEV | K | K | K |

K | K | LPL | lysophospholipids |

NADPH | reduced form of nicotinamide adenine dinucleotide phosphate | NIA1 | nitrate reductase |

NO | nitric oxide | OST1 | protein kinase open stomata 1 |

PA | phosphatidic acid | PEPC | phosphoenolpyruvate carboxylase |

phot1 | phototropin 1 | phot1 | 14-3-3 protein bound phototropin 1 |

phot2 | phototropin 2 | Photophos-phorylation | light-dependent reactions of photosynthesis |

PIP2 | phosphatidylinositol 4,5-bisphosphate located in the cytosol | PIP2 | phosphatidylinositol 4,5-bisphosphate located at the plasma membrane |

PLA | phospholipase A2β | PLC | phospholipase C |

PLD | phospholipase D | PMV | electric potential difference across the plasma membrane |

PP1 | the catalytic subunit of type 1 phosphatase located in the nucleus | PP1 | the catalytic subunit of type 1 phosphatase located in the cytosol |

protein kinase | a serine/threonine protein kinase that directly phosphorylates the plasma membrane H-ATPase | PRSL1 | type 1 protein phosphatase regulatory subunit 2-like protein1 |

RIC7 | ROP-interactive CRIB motif-containing protein 7 | ROP2 | small GTPase ROP2 |

ROS | reactive oxygen species | [Ca | cytosolic Ca |

[Cl | cytosolic/vacuolar Cl | [K | cytosolic/vacuolar K |

[malate | apoplastic/ cytosolic/vacuolar malate | [NO | apoplastic/cytosolic/vacuolar nitrate concentration |

## Results

### Network reduction

The Sun et al. model has a huge state space of ~10^{31} states, making its analysis difficult. To obtain a smaller state space, we reduce the size of the network by applying a network reduction technique developed by Saadatpour et al. [23] that is proven to preserve the attractors of Boolean models (see Methods). All source nodes other than the four signals (blue light, red light, CO_{2}, and abscisic acid) and all simple mediator nodes are identified and reduced. This process is done iteratively until it cannot be done any more. A total of 7 source nodes (14-3-3 protein_{phot1}, PIP2_{C}, AtNOA1, Nitrate, PP1_{cn}, mitochondria, and CHL1), and 19 simple mediator nodes (phot1, phot2, NIA1, H^{+}-ATPase, LPL, ATP, acid. of apoplast, [NO_{3}
^{−}]_{v}, [Cl^{−}]_{v}, NADPH, [malate^{2−}]_{v}, PA, ABA receptors, OST1, PRSL1, PIP2_{PM}, AtrbohD/F, Nitrite, and phot1_{complex}) are eliminated. Several of the simple mediator nodes form linear paths (e.g. phot1, OST1) thus their iterative reduction shortens the linear paths in the network. In addition, 16 of the 19 reduced mediators have a regulatory function of the form “B* = A”. It is intuitive that reduction of this node type preserves the attractors.

We do not eliminate the four signal nodes because we want to simultaneously explore all the combinations of input signals. We also choose to not reduce the five nodes (K_{in}, K_{out}, K_{c}, Ca^{2+}-ATPase, mesophyll cell photosynthesis) whose merging with their sole regulator would result in a self-loop (self-regulation), because such self-loops may be difficult to interpret. Two additional nodes with significant biological meaning to the network (sucrose, stomatal opening), are not reduced either.

Another form of network reduction is the elimination of redundant edges (see Methods). After removal of redundant edges, the node CDPK becomes a sink node, thus it can also be eliminated. The reduction of the above-described nodes and redundant edges simplifies the network from 70 nodes to 42 nodes, with an estimated state space of ~10^{22} states.

### Simplification of regulatory functions

In order to further reduce the state space from ~10^{22} to a manageable size, we grouped state values so that nodes are represented with fewer states. This grouping was guided by the 66 experimental observations summarized in Sun et al.; we aimed to maintain the reduced model’s results consistent with these experimental observations.

^{+}]

_{v}due to charge balance. The anion contributions are [malate

^{2−}]

_{v}contribution ≤ 0.425 × [K

^{+}]

_{v}; [NO

_{3}

^{−}]

_{v}contribution ≤0.10 × [K

^{+}]

_{v}; [Cl

^{−}]

_{v}contribution ≤ 0.05 × [K

^{+}]

_{v}. The primary contributions come from [K

^{+}]

_{v}and sucrose. We grouped the stomatal opening values into 6 groups with different [K

^{+}]

_{v}and sucrose values (see Table 2 and Additional file 1).

Grouping of the stomatal opening values by the level of [K^{+}]_{v} and sucrose

[K | Sucrose | Stomatal opening value in the Sun et al. model | Simplified stomatal opening value |
---|---|---|---|

0 | 0 | 0 | 0 |

0 | 1 or 2 | 1 or 2 | 1 |

1 | 0 | 1.58 | 1 |

1.8 | 1 | 3.84 | 2 |

1.5 | 2 | 4.36 | 2 |

2 | 0 or 1 | 3.15 or 4.15 | 3 |

4.5 | 0 or 2 | 5.18 or 8.92 | 3 |

6 | 0 | 9.28 or 9.45 | 5 |

6 | 2 | 11.28 or 11.45 | 5 |

9 | 0 or 2 | 14.01 or 16.01 | 6 |

The first two columns indicate the [K^{+}]_{v} and sucrose levels. The third column is the possible values of stomatal opening in the Sun et al. model for the given [K^{+}]_{v} and sucrose levels. Note that here we only show [K^{+}]_{v}, sucrose and stomatal opening value combinations observed in the simulations of the 66 experimentally studied scenarios reported by Sun et al. [15]. More stomatal opening values are possible when considering node perturbations. The 4th column shows the simplified stomatal opening level after grouping. The update function for the simplified stomatal opening level covers all possible values of [K^{+}]_{v} and sucrose (see Additional file 1).

Similarly to the original model, the simplified states represent qualitative, relative categories. For example, a stomatal opening level of 2 is not twice as high as level 1. We choose the simplified stomatal opening values so that there is no state “4”, to better reflect an experimentally observed synergistic effect between blue and red light [18, 19, 27]. Simulation results with the simplified regulatory function are that under monochromatic red light stomatal opening =1; under monochromatic blue light stomatal opening =3; under dual beam the stomatal opening =5, which is larger than the sum “1 + 3”. This qualitatively reproduces the experimental observation that under dual beam illumination stomata open to a size much larger than the sum of opening under monochromatic blue or red light.

We find by simulation of the reduced model, using the same initial condition as the Sun et al. model, that the simplification of the stomatal opening regulatory function results in only 3 additional cases of inconsistency with experimental observations out of a total of 66 experimentally studied scenarios. Additional file 2 lists all experimental observations and compares them to the relevant simulation results. Ignoring the contribution of malate^{2−}, NO_{3}
^{−}, and RIC7 to stomatal opening each causes one additional discrepancy; ignoring Cl^{−} does not cause any additional discrepancy. Ignoring these nodes trades a decrease in accuracy for a significant increase in simplicity.

^{2−}]

_{a}, [malate

^{2−}]

_{c}, starch, [Cl

^{−}]

_{c}, [NO

_{3}

^{−}]

_{c}, [NO

_{3}

^{−}]

_{a}, ROP2, RIC7, ABC, and PEPC. The only edge from these nodes to other nodes is [malate

^{2−}]

_{a}→ AnionCh. In section 3 of Additional file 3 we show that eliminating this edge does not change the system’s long-term behavior, i.e. attractors. Also, the regulatory function describing the cytosolic K

^{+}concentration, [K

^{+}]

_{c}, can be simplified without loss, as described in section 3 of Additional file 3. After this simplification we have a network of 32 nodes, 81 edges, indicated on Fig. 2. We will refer to this model as the “reduced model”. A list of nodes and their regulatory functions is provided in Additional file 1.

Identifying strongly connected components (SCCs) is important for attractor analysis, as complex dynamic behavior such as oscillations or multi-stability requires feedback loops [7]. There are three SCCs in the network of the reduced model, as marked in Fig. 2. The NO cycle contains three nodes and three positive edges. The C_{i} SCC contains three nodes, which form two negative feedback loops. The Ion SCC is the most complex, containing 13 nodes and 26 edges, 7 of which are negative.

Next we perform attractor analysis using two methods: 1. by converting the reduced model to Boolean and applying two analysis tools; 2. by analyzing the regulatory functions theoretically. The former method finds all stable steady states and candidate oscillations; the latter confirms the results of the first method and gives insight about perturbation scenarios.

### Conversion of nodes from multi-level to Boolean states and attractor analysis

Example of Boolean conversion

Level of the original node | State of Boolean node_2 | State of Boolean node_1 |
---|---|---|

0 | 0 | 0 |

1 | 0 | 1 |

2 | 1 | 0 |

3 | 1 | 1 |

More detailed examples of the conversion of the states and regulatory function of specific nodes are given in the Additional file 4. We will refer to the reduced model after conversion to Boolean variables as the “Boolean-converted reduced model”. The regulatory functions of the Boolean-converted reduced model are available in Additional file 5. When simulating the Boolean-converted reduced model, all the Boolean nodes that represent the same entity (the same multi-level node) are updated simultaneously. In this way the state transitions of the reduced model will be kept the same in the Boolean-converted reduced model, and therefore the Boolean conversion will not cause additional discrepancies from experimental observations.

_{2}, CO

_{2}_high). We find two possible stable motifs, corresponding to the self-regulatory node PMV_pos (one of the two Boolean nodes associated with the multi-level node PMV, see Additional files 4 and 5), in conditions where the H

^{+}-ATPase

_{complex}is inactive. These two stable motifs indicate the bistability of PMV. Under its influence, another node, K

_{out}, will also be bistable. The algorithm also indicates that for any signal combination, every node, except [Ca

^{2+}]

_{c}and Ca

^{2+}-ATPase, will stabilize in a fixed state. [Ca

^{2+}]

_{c}has three states, and in the Boolean-converted model it is represented by two nodes, Cac and Cac_high. Cac_high, which represents the higher level of [Ca

^{2+}]

_{c}, stabilizes at zero in all situations. Cac and Ca

^{2+}-ATPase may oscillate in conditions where blue light is present and ABA is absent (a total of six cases, two of which allow PMV bistability). Table 4 summarizes key features of the attractors found by the stable motif algorithm for all 24 input combinations. Attractors where Ca

^{2+}oscillation is not possible are fixed points (stable steady states).

Summary of the attractors found using the stable motif algorithm

BL | RL | CO | CO | ABA | SO (Bool) | SO | Ca | PMV_pos bistability |
---|---|---|---|---|---|---|---|---|

0 | 0 | Any | Any | Any | 000 | 0 | No | Yes |

0 | 1 | 0 | 0 | 1 | 000 | 0 | No | No |

0 | 1 | 1 | Any | 1 | 000 | 0 | No | Yes |

1 | Any | 1 | 0 | 1 | 000 | 0 | No | No |

1 | Any | 1 | 1 | 1 | 000 | 0 | No | Yes |

0 | 1 | 1 | Any | 0 | 010 | 1 | No | Yes |

1 | Any | 1 | 1 | 0 | 010 | 1 | Yes | Yes |

0 | 1 | 0 | 0 | 0 | 101 | 3 | No | No |

1 | 0 | 1 | 0 | 0 | 101 | 3 | Yes | No |

1 | Any | 0 | 0 | 1 | 101 | 3 | No | No |

1 | 0 | 0 | 0 | 0 | 110 | 5 | Yes | No |

1 | 1 | 1 | 0 | 0 | 110 | 5 | Yes | No |

1 | 1 | 0 | 0 | 0 | 111 | 6 | Yes | No |

We verified the obtained attractors with GINsim [12], a software suite capable of model construction, simulation, and analysis. GINsim can compute all stable steady states (called stable states in GINsim), or determine complex attractors by mapping the state transitions. The stable steady states found by GINsim are identical to those found by the stable motif algorithm. To verify and further explore the complex attractors, we use the simulation function of GINsim, starting from a state in the complex attractor. The result that the system oscillates between four states, where only the state of Cac and Ca^{2+}-ATPase changes, agrees with the findings of the stable motif algorithm. We summarize the GINsim computation/simulation results in Additional file 7. Additional file 8 indicates the Boolean-converted reduced model in SBML-qual format [29], a general format for biological model to be analyzed using various tools including GINsim.

We can also connect the stable motif analysis results to network reduction. We have previously decided to not reduce the four nodes that correspond to input signals. If we do consider a specific input combination when using network reduction, e.g. blue light and red light with normal CO_{2} without ABA, we can reduce much more of the network: two of the three SCCs, namely the NO cycle and the C_{i} SCC, will stabilize and can be eliminated. Only the Ion SCC and its sole output stomatal opening remain, indicating that this SCC is not driven solely by the external signals and has the capacity for oscillations or multi-stability. This is consistent with the results found by stable motif analysis, according to which the NO cycle and the C_{i} SCC attain a steady state and the Ion SCC admits a [Ca^{2+}]_{c} - Ca^{2+}-ATPase oscillation and PMV bistability. This consistency supports the appropriateness of the network reduction method and of the Boolean conversion.

### Theoretical analysis of the reduced model

To gain additional insight into the attractors of the reduced model and their potential changes due to node perturbations, we analyze the reduced model theoretically. Specifically, we aim to answer the question: Can there be other types of oscillation, or can there be additional multi-stability, if a node is knocked out (fixed in the OFF state) or is constitutive active (fixed in the highest state)?

We first test whether the network and regulatory rules allow multi-stability or oscillations. This analysis is based on R. Thomas’s conjectures [7]: The presence of a positive (negative) feedback loop - a cycle with an even (odd) number of inhibitory edges - in the network is a necessary but not sufficient condition for the occurrence of multiple steady states (oscillations). The conjectures have been proven in the case of discrete dynamic systems [31–34]. Since only feedback loops are candidates for potential multi-stability or oscillations, we analyze the regulatory functions of each strongly connected component of the network. For each feedback loop, we identify a sufficient condition for the nodes to stabilize in a specific state. The violation of this condition becomes a further necessary condition of multi-stability or oscillation. Here we describe the main steps and results of the analysis; the detailed analysis is in Additional file 3.

The NO cycle is composed of the nodes PLD, ROS, NO, and the three positive edges between them. It does not have any negative edges, so it cannot oscillate. A fixed ABA value is sufficient to stabilize each node of the cycle in a specific state, thus the cycle does not admit multi-stability under any perturbation.

The C_{i} SCC has three nodes, C_{i}, mesophyll cell photosynthesis (MCPS), carbon fixation, and four edges that form two negative feedback loops, one between carbon fixation and C_{i}, and the other between C_{i} and MCPS. Despite the existence of negative feedback, this cycle will stabilize if given a fixed CO_{2} value. From this we know that this cycle cannot oscillate or admit multi-stability under any perturbation.

^{2+}]

_{c}, which has states 0,1, and 2, cannot enter state 2 in the long term under any perturbation. Since most nodes respond to [Ca

^{2+}]

_{c}only if [Ca

^{2+}]

_{c}=2, we can eliminate all edges that depend only on “[Ca

^{2+}]

_{c}=2”, and obtain a simplified Ion SCC, as shown in Fig. 3. The Ca

^{2+}SCC ([Ca

^{2+}]

_{c}, Ca

^{2+}ATPase, PLC, CaR) now becomes a sink SCC. The only negative edge in this sub-network is from Ca

^{2+}-ATPase to [Ca

^{2+}]

_{c}. These two nodes are known to oscillate. The positive feedback loop formed by [Ca

^{2+}]

_{c}, PLC, and CaR will stabilize if given fixed inputs. So there cannot be multi-stability. For the nodes outside of the Ca

^{2+}feedback loops, we show that the edges from KEV and [K

^{+}]

_{v}are redundant in the long term, so there are no feedback loops except the PMV self-loop. PMV is not capable of having oscillations, but can have bistability (as also indicated by the stable motif analysis). The bistability can affect at most one other node, K

_{out}, under any perturbation. This means that the bistability has very limited effect on the attractor of the reduced model.

Now we can summarize our conclusions and return to the question we sought to answer: there is no oscillation except in the calcium nodes; there is no multi-stability except in the nodes PMV and K_{out}. These statements are true under any perturbation. Moreover, for the calcium oscillation, [Ca^{2+}]_{c} cannot enter the state 2, so the sub-network between [Ca^{2+}]_{c} and Ca^{2+}-ATPase is a negative feedback loop between two Boolean nodes, with the regulatory functions Ca^{2+} ATPase^{*} = [Ca^{2+}]_{c}; [Ca^{2+}]_{c}* = not Ca^{2+} ATPase. It results in the simplest type of oscillation, as also found by GINsim simulation. For the PMV bistability, even if the bistability exists, most nodes, especially the output node stomatal opening, still have a unique value. Thus the theoretical analysis, in agreement with the computational analysis, leads to very strong conclusions about the reduced model’s dynamic repertoire.

We can also show that the reduction or Boolean conversion did not change the attractors of the Sun et al. model. Although the reduction we used is only proven in the Boolean case, Naldi et al. showed that for multi-valued models, removal of non-autoregulated nodes, like in our reduction, preserves crucial dynamical properties [35], including fixed point attractors and the two-node simple oscillation we found. So our reduction is valid in this specific model. To confirm that the Boolean conversion preserved attractors, we note that in the Boolean-converted reduced model we found fixed point attractors and a complex attractor in which only two nodes oscillate. Because the only potential change to attractors as a consequence of the conversion is merging of complex attractors [26], it is straightforward that the attractors have been conserved during the conversion, as the two-node oscillation found is the simplest type of complex attractor and cannot be a result of attractor merging. In addition, using general asynchronous update instead of random order asynchronous update does not cause any changes to the attractor, because the update schemes do not affect fixed points or the two-node simple oscillation we found.

### Stability of guard cell signal transduction

_{2}and ABA conditions. For each signal combination, we set the perturbed node’s initial state and regulatory function to 0, initialize the rest of the nodes in the condition representative of closed stomata, and then simulate the reduced model until it reaches its attractor. In the absence of ABA under each light and CO

_{2}condition, 60–90 % perturbation scenarios produce the same stomatal opening value as the unperturbed system (Table 5). These results are similar to those reported by Sun et al. for the original model [15] (see Additional file 9). In the presence of ABA 50–90 % perturbation scenarios produce the same stomatal opening value as the unperturbed system, and 4–16 % knockouts lead to a higher stomatal opening value. Perturbations in the ABA = 1 case were not studied by Sun et al., but our simulations of the original model give the same qualitative results as the reduced model. These results indicate the closeness of the perturbed attractor (at least in terms of the stomatal opening value) to the unperturbed attractor in more than 50 % of single node perturbations. They also suggest the resilience of the stomatal opening process against internal failures and perturbations.

Summary of systematic perturbation results

Light, CO | Unperturbed SO level | Simplified SO level | Percentage of cases with unchanged SO value | |||||||
---|---|---|---|---|---|---|---|---|---|---|

0 | 1 | 2 | 3 | 5 | 6 | |||||

Percentage of single knockouts that lead to each SO level | ||||||||||

Dual Beam | Mod. CO | ABA OFF | 5 | 4 % | 31 % | 65 % | 65 % | |||

Low CO | 6 | 31 % | 4 % | 65 % | 65 % | |||||

High CO | 1 | 4 % | 96 % | 96 % | ||||||

Blue Light | Mod. CO | 3 | 35 % | 65 % | 65 % | |||||

Low CO | 5 | 31 % | 4 % | 65 % | 65 % | |||||

High CO | 1 | 4 % | 96 % | 96 % | ||||||

Red Light | Mod. CO | 1 | 4 % | 96 % | 96 % | |||||

Low CO | 3 | 35 % | 65 % | 65 % | ||||||

High CO | 1 | 4 % | 96 % | 96 % | ||||||

Dual Beam | Mod. CO | ABA ON | 0 | 85 % | 4 % | 8 % | 4 % | 85 % | ||

Low CO | 3 | 46 % | 50 % | 4 % | 50 % | |||||

Blue Light | Mod. CO | 0 | 85 % | 4 % | 8 % | 4 % | 85 % | |||

Low CO | 3 | 46 % | 50 % | 4 % | 50 % | |||||

Red Light | Low CO | 0 | 96 % | 4 % | 96 % |

### Extending the conclusions to the original model

We found that in the reduced model there is no oscillation except in the calcium nodes; there is no multi-stability except in the nodes PMV and K_{out}. Because the reduction we used has been shown to conserve attractors [23, 35], we know that our attractor conclusions can be immediately extended to all nodes in the original model except the reduced nodes and stomatal opening. Next we extend the attractor analysis to include the reduced nodes as well.

First we consider the nodes reduced during the first step of network reduction, i.e. non-signal source nodes and simple mediator nodes. These nodes are trivially incapable of having multi-stability and oscillations themselves, so we need only to consider their perturbations. Perturbation of a simple mediator node can always be replaced by a corresponding (set of) perturbation(s) in the mediator node’s direct successor(s), so these perturbations have already been considered. Perturbing a non-signal source node may theoretically cause a difference, however the nodes in this category in the Sun et al. model represent molecules that are abundant in the cell or cell environment, thus their perturbation is not biologically relevant or practical.

Next we consider the anion nodes reduced due to the simplified stomatal opening rule. Recall that these nodes do not affect other nodes except stomatal opening in the long term. There cannot be multi-stability in anion nodes unless the assumptions of sufficient initial [NO_{3}
^{−}]_{a} and starch concentration, and sufficient initial mitochondrial TCA cycle activity are violated (details are provided in Additional file 3, section 5 and 6). Since there is no support for interventions that would lead to the violation of these assumptions, it is reasonable to conclude that no multi-stability can be found in the reduced nodes under biologically relevant situations. We also found that there can be an additional oscillation in the RIC7 path (involving the nodes ROP2, RIC7 and SO) when a special set of perturbations is applied. Under that case, the nodes RIC7 and SO will oscillate. Since the effect of this behavior is small (within 5 % of the unperturbed SO value in the Sun et al. model [15]), it has little biological significance. There are no more possible oscillations as there are no more negative feedback loops. To conclude, the original Sun et al. model has oscillations only in cytosolic Ca^{2+} ([Ca^{2+}]_{c}) and Ca^{2+} ATPase, and has multi-stability only in PMV and K_{out}, under situations that are biologically meaningful.

## Discussion

The conclusions we obtained can tell us how to control this network model. Generally in engineering applications, control means to drive a system into an arbitrary state [36, 37]. However in biological systems, it is more meaningful to drive the system into one of its natural attractors rather than into an arbitrary state, as the attractors correspond to stable phenotypes [38]. To control the attractor of a Boolean system, one needs to control only its input nodes and a subset of nodes in each stable motif [39]. Our integrated analysis, involving Boolean conversion, indicates that to control the attractor that the stomatal opening network evolves into, one only needs to control the input signals and PMV, even in case of perturbations. In particular, to control the stomatal opening value, one only needs to control the input signals, under any perturbation.

Nodes whose knockouts diminish ABA’s inhibition of stomatal opening

Light, CO | Unperturbed SO level | Nodes whose knockout results in a partially restored SO, and the corresponding SO value | |||||
---|---|---|---|---|---|---|---|

CO | NO | PLD | ROS | AnionCh | |||

Dual Beam | Moderate CO | 0 | 3 | 3 | 5 | 3 | 2 |

Blue Light | 0 | 3 | 2 | 3 | 2 | 1 | |

Red Light | 0 | 3 | 1 |

Our results reproduce the observation that knockout of nodes in the ABA pathway (PLD, NO, ROS) can cause partial reversals of ABA’s effect. We find that AnionCh knockout can partially restore stomatal opening inhibited by ABA, a result not reported by Sun et al., but which is supported by experimental evidence [40]. In addition, Table 6 offers a new biological prediction: low CO_{2} concentration can partially restore stomatal opening when ABA is present. This is consistent with the knowledge that CO_{2}-free air promotes stomatal opening in the absence of ABA [41]. This CO_{2} effect suggests a mechanism of cross-talk between CO_{2} and ABA. Importantly, apart from the five nodes listed in Table 6, no other node’s knockout can reverse ABA’s inhibition of stomatal opening. The perturbation results of Table 5 offer many more new predictions.

Our combination of techniques offers a powerful framework for determining the dynamic repertoire of a multi-level dynamic model. Multi-level models are more accurate than Boolean models in describing the quantitative characteristics of dynamic systems, but there are few general methods to analyze multi-level models [10, 12]. By combining different existing methods, we were able to overcome the limitations of each method. Our successful combination of existing methods offers a promising way to analyze multi-level models, and might point towards a general strategy to analyze the attractors of multi-level models, biological or non-biological.

A notable future direction for this work is to develop an alternative way to determine the attractors of multi-level models by extending the concept of stable motifs. Compared with conversion to a Boolean model, then applying Boolean stable motif algorithm, extending the stable motif algorithm to multi-level models can avoid potential attractor change issues. Development of such a technique will allow easy and powerful attractor analysis for multi-level models.

## Conclusions

We obtained a very strong conclusion about the attractors of the Sun et al. stomatal opening model: under any combination of sustained signals, all nodes in the model converge into steady states, with the potential exception of the cytosolic Ca^{2+} ([Ca^{2+}]_{c}) and Ca^{2+} ATPase. Variations in the initial condition of non-source nodes or in process timing (node update sequence) can drive at most two nodes, PMV and K_{out}, into a different attractor. This high degree of attractor similarity is somewhat unexpected, as the network has a large strongly connected component and several feedback loops. Thus, despite the decidedly non-linear structure of the network, most parts of the system behave in the consistent manner of a linear pathway. This is a distinct feature of the stomatal opening model: many dynamic models of biological systems have multiple, diverse attractors [22, 42]. The models of these systems will evolve into drastically different attractors when starting from different initial conditions, sometimes even when starting from the same initial condition, demonstrating different biological trajectories. In the stomatal opening model, however, the uniqueness of the steady state stomatal opening level suggests that the final extent of the stomatal opening response is robust and resilient against changes in initial conditions or in timing. Note that although a change in the initial condition will not change the steady-state opening level, it may change the steady state of PMV and K_{out}, and may change how fast the system converges to an attractor.

We also showed that the reduced stomatal opening model does not admit additional, emergent oscillations or multi-stability under any biologically relevant node perturbation (knockout or constitutive activity). We further demonstrate the robustness of the system by examining the stomatal opening level under single node knockouts: in most cases the signals are still likely to propagate and lead to a similar degree of stomatal opening as in the absence of perturbation. This robustness is unlike a single linear pathway, which would be very sensitive to node disruption. We suggest that the role of the strongly connected components in the network could be to provide multiple paths for the signal to propagate, but at the same time not allowing extensive multistability or oscillations. Our innovative combination of existing methods offers a promising way to analyze multi-level models.

## Declarations

### Acknowledgements

The authors thank Zhongyao Sun, Jorge G. T. Zañudo and Prof. Sarah Assmann for helpful discussions.

### Funding

This project was supported by National Science Foundation (NSF) grants IIS 1160995, PHY 1205840 and MCB 1244303. The NSF had no role in the design of the study, analysis and interpretation of data or in writing the manuscript.

### Availability of data and materials

The datasets supporting the conclusions of this article are included within the article and its additional files.

### Authors’ contributions

XG developed the reduced model and performed the analysis under the advice and supervision of RA. Both authors wrote the manuscript. Both authors have read and approved the final version of the manuscript.

### Competing interests

The authors declare that they have no competing interests.

### Consent for publication

Not applicable. The article does not contain any individual person’s data.

### Ethics approval and consent to participate

Not applicable. The article does not involve human participants, human data or human tissue.

**Open Access**This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

## Authors’ Affiliations

## References

- Stigler B, Chamberlin HM. A regulatory network modeled from wild-type gene expression data guides functional predictions in Caenorhabditis elegans development. BMC Syst Biol. 2012;6:77.View ArticlePubMedPubMed CentralGoogle Scholar
- Chifman J, et al. The core control system of intracellular iron homeostasis: a mathematical model. J Theor Biol. 2012;300:91–9.View ArticlePubMedPubMed CentralGoogle Scholar
- Massague J. TGF-beta signal transduction. Annu Rev Biochem. 1998;67:753–91.View ArticlePubMedGoogle Scholar
- Xu HL, et al. Construction and Validation of a Regulatory Network for Pluripotency and Self-Renewal of Mouse Embryonic Stem Cells. PLoS Comput Biol. 2014;10(8):e1003777.Google Scholar
- Kestler HA, et al. Network modeling of signal transduction: establishing the global view. Bioessays. 2008;30(11–12):1110–25.View ArticlePubMedGoogle Scholar
- Tyson JJ, Chen K, Novak B. Network dynamics and cell physiology. Nat Rev Mol Cell Biol. 2001;2(12):908–16.View ArticlePubMedGoogle Scholar
- Thomas R, European Molecular Biology Organization. Kinetic logic : a Boolean approach to the analysis of complex regulatory systems : proceedings of the EMBO course "Formal analysis of genetic regulation," held in Brussels, September 6–16, 1977, Lecture notes in biomathematics. Berlin; New York: Springer; 1979. p. xiii, 507.Google Scholar
- Kauffman SA. Metabolic stability and epigenesis in randomly constructed genetic nets. J Theor Biol. 1969;22(3):437.View ArticlePubMedGoogle Scholar
- Glass L, Kauffman SA. Logical analysis of continuous, nonlinear biochemical control networks. J Theor Biol. 1973;39(1):103–29.View ArticlePubMedGoogle Scholar
- Wang RS, Saadatpour A, Albert R. Boolean modeling in systems biology: an overview of methodology and applications. Phys Biol. 2012;9(5):055001.View ArticlePubMedGoogle Scholar
- Miskov-Zivanov N, et al. The duration of T cell stimulation is a critical determinant of cell fate and plasticity. Sci Signal. 2013;6(300):ra97.View ArticlePubMedPubMed CentralGoogle Scholar
- Chaouiya C, Naldi A, Thieffry D. Logical modelling of gene regulatory networks with GINsim. Methods Mol Biol. 2012;804:463–79.View ArticlePubMedGoogle Scholar
- Deritei D, et al. Principles of dynamical modularity in biological regulatory networks. Sci Rep. 2016;6:21957.View ArticlePubMedPubMed CentralGoogle Scholar
- Murrugarra D, Laubenbacher R. Regulatory patterns in molecular interaction networks. J Theor Biol. 2011;288:66–72.View ArticlePubMedGoogle Scholar
- Sun Z, et al. Multi-level modeling of light-induced stomatal opening offers new insights into its regulation by drought. PLoS Comput Biol. 2014;10(11):e1003930.View ArticlePubMedPubMed CentralGoogle Scholar
- Li S, Assmann SM, Albert R. Predicting essential components of signal transduction networks: a dynamic model of guard cell abscisic acid signaling. PLoS Biol. 2006;4(10):e312.View ArticlePubMedPubMed CentralGoogle Scholar
- Schroeder JI, et al. Guard Cell Signal Transduction. Annu Rev Plant Physiol Plant Mol Biol. 2001;52:627–58.View ArticlePubMedGoogle Scholar
- Shimazaki K, et al. Light regulation of stomatal movement. Annu Rev Plant Biol. 2007;58:219–47.View ArticlePubMedGoogle Scholar
- Assmann SM. Enhancement of the Stomatal Response to Blue Light by Red Light, Reduced Intercellular Concentrations of CO(2), and Low Vapor Pressure Differences. Plant Physiol. 1988;87(1):226–31.View ArticlePubMedPubMed CentralGoogle Scholar
- Bergmann DC, Sack FD. Stomatal development. Annu Rev Plant Biol. 2007;58:163–81.View ArticlePubMedGoogle Scholar
- MacArthur BD, Ma'ayan A, Lemischka IR. Systems biology of stem cell fate and cellular reprogramming. Nat Rev Mol Cell Biol. 2009;10(10):672–81.PubMedPubMed CentralGoogle Scholar
- Steinway SN, et al. Network modeling of TGFbeta signaling in hepatocellular carcinoma epithelial-to-mesenchymal transition reveals joint sonic hedgehog and Wnt pathway activation. Cancer Res. 2014;74(21):5963–77.View ArticlePubMedPubMed CentralGoogle Scholar
- Saadatpour A, Albert R, Reluga TC. A Reduction Method for Boolean Network Models Proven to Conserve Attractors. SIAM J Appl Dyn Syst. 2013;12(4):1997–2011.View ArticleGoogle Scholar
- Ansotegui C, Manya F. Mapping problems with finite-domain variables to problems with Boolean variables. Theory Appl of Satisfiability Test. 2005;3542:1–15.View ArticleGoogle Scholar
- Van Ham P. How to deal with more than two levels. In: Thomas R, editor. Kinetic logic : a Boolean approach to the analysis of complex regulatory systems : proceedings of the EMBO course "Formal analysis of genetic regulation," held in Brussels, September 6–16, 1977. Berlin; New York: Springer; 1979. p. 326–44.Google Scholar
- Didier G, Remy E, Chaouiya C. Mapping multivalued onto Boolean dynamics. J Theor Biol. 2011;270(1):177–84.View ArticlePubMedGoogle Scholar
- Karlsson PE. Blue light regulation of stomata in wheat seedlings. I. Influence of red background illumination and initial conductance level. Physiol Plant. 1986;66:5.Google Scholar
- Zanudo JG, Albert R. An effective network reduction approach to find the dynamical repertoire of discrete dynamic networks. Chaos. 2013;23(2):025111.View ArticlePubMedGoogle Scholar
- Chaouiya C, et al. SBML qualitative models: a model representation format and infrastructure to foster interactions between qualitative modelling formalisms and tools. BMC Syst Biol. 2013;7:135.View ArticlePubMedPubMed CentralGoogle Scholar
- Veliz-Cuba A, et al. Steady state analysis of Boolean molecular network models via model reduction and computational algebra. BMC Bioinformatics. 2014;15:221.View ArticlePubMedPubMed CentralGoogle Scholar
- Remy E, Ruet P, Thieffry D. Graphic requirements for multistability and attractive cycles in a Boolean dynamical framework. Adv Appl Math. 2008;41(3):335–50.View ArticleGoogle Scholar
- Remy E, Ruet P. On differentiation and homeostatic behaviours of Boolean dynamical systems. Lect Notes Bioinformatics. 2007;4780:92–101.Google Scholar
- Richard A, Comet J-P. Necessary conditions for multistationarity in discrete dynamical systems. Discret Appl Math. 2007;155(18):2403–13.View ArticleGoogle Scholar
- Richard A. Negative circuits and sustained oscillations in asynchronous automata networks. Adv Appl Math. 2010;44(4):378–92.View ArticleGoogle Scholar
- Naldi A, et al. Dynamically consistent reduction of logical regulatory graphs. Theor Comput Sci. 2011;412(21):2207–18.View ArticleGoogle Scholar
- Liu YY, Slotine JJ, Barabasi AL. Controllability of complex networks. Nature. 2011;473(7346):167–73.View ArticlePubMedGoogle Scholar
- Lin CT. Structural controllability. IEEE Trans Autom Control. 1974;AC19(3):201–8.Google Scholar
- Mochizuki A, et al. Dynamics and control at feedback vertex sets. II: A faithful monitor to determine the diversity of molecular activities in regulatory networks. J Theor Biol. 2013;335:130–46.View ArticlePubMedGoogle Scholar
- Zanudo JG, Albert R. Cell fate reprogramming by control of intracellular network dynamics. PLoS Comput Biol. 2015;11(4):e1004193.View ArticlePubMedPubMed CentralGoogle Scholar
- Schwartz A, et al. Anion-Channel Blockers Inhibit S-Type Anion Channels and Abscisic Acid Responses in Guard Cells. Plant Physiol. 1995;109(2):651–8.PubMedPubMed CentralGoogle Scholar
- Kim TH, et al. Guard cell signal transduction network: advances in understanding abscisic acid, CO2, and Ca2+ signaling. Annu Rev Plant Biol. 2010;61:561–91.View ArticlePubMedPubMed CentralGoogle Scholar
- Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J Theor Biol. 2003;223(1):1–18.View ArticlePubMedGoogle Scholar