The spruce budworm and forest: a qualitative comparison of ODE and Boolean models

Boolean and polynomial models of biological systems have emerged recently as viable companions to differential equations models. It is not immediately clear however whether such models are capable of capturing the multi-stable behaviour of certain biological systems: this behaviour is often sensitive to changes in the values of the model parameters, while Boolean and polynomial models are qualitative in nature. In the past few years, Boolean models of gene regulatory systems have been shown to capture multi-stability at the molecular level, confirming that such models can be used to obtain information about the system’s qualitative dynamics when precise information regarding its parameters may not be available. In this paper, we examine Boolean approximations of a classical ODE model of budworm outbreaks in a forest and show that these models exhibit a qualitative behaviour consistent with that derived from the ODE models. In particular, we demonstrate that these models can capture the bistable nature of insect population outbreaks, thus showing that Boolean models can be successfully utilized beyond the molecular level.


Introduction
The spruce budworm, Choristoneura fumiferana (Clements), is the most serious defoliator of spruce (Picea spp.) and balsam fir (Abies balsamea) trees in the boreal forest of North America. The species is transcontinental, occurring from Newfoundland to Alaska and down into the northeastern US (Lyons et al., 1997). Due to its high impact on the forest resources, much research has been devoted to developing different types of models and simulators for the budworm-forest system (e.g. James et al., 2010;Ludwig et al., 1978;Lyons et al., 1997;Régnière & You, 1991) utilizing both stochastic and deterministic approaches, as well as tools for screening for potential outbreaks and methods for control (see, e.g. Nealis & Nault, 2005). The outbreak dynamics and consequences (ecological and economic) have also been linked to climate change (Jamieson et al., 2012;Thomson, 2009).
It is well known that the budworm population could be either very low or very high under similar external conditions (Morris, 1963). This property of a system is called bistability.
A bistable system can reside in one of the two states for a long time unless a strong external perturbation forces a change. Bistability also provides a discontinuous switching among the stable steady states, possibly converting graded inputs into switch-like outputs. Bistability comes with hysteresis, which can be defined as dependence of a steady state of a system not only on its current state but also on its past state(s). To predict a future state of a bistable system, both its current state and its history must be known. This notion goes back to the classical work by Landau (1937) in which he develops a qualitative theory of phase transitions.
In this paper, we study Boolean approximations of a classical ordinary differential equations (ODE) model of the budworm-forest dynamics that was developed and analysed in Ludwig et al. (1978). The ODE model under consideration captures the bistability of the system and our goal was to show that Boolean models can also capture that feature. That Boolean models can be bistable was first established in Veliz- Cuba and Stigler (2011), followed by work on Boolean approximations of ODE models that preserve bistability (Hinkelmann & Laubenbacher, 2011;Robeva & Yildirim, 2013). All of these results however were obtained for systems at the molecular level where much of the system's dynamics are known a priori to be close to Boolean (e.g. genes are either on or off, proteins are either made or not, etc.). Thus, a question regarding the general suitability of Boolean models at the population level deserves attention and the ability of such models to capture biologically meaningful states of bistability, as in the case of the budworm-forest system, is of particular interest. Answering these questions was the main motivator for the work presented here. For the Boolean models, we have purposefully used qualitative analyses that parallel those utilized to study the ODE models in Ludwig et al. (1978). In our opinion, this approach emphasizes the high level of similarity with which qualitative explorations of these different types of models can be carried out.
The paper is organized as follows. In Section 2, we present in some detail the ODE models introduced in Ludwig et al. (1978) together with their qualitative analyses. In Section 3, we present our Boolean approximations of these models and show that they exhibit essentially the same qualitative behaviour, including bistability. In Section 4, we discuss possible generalizations of the Boolean models and the robustness of their outcomes, comparing with results obtained under different time-update schedules.

ODE models of spruce budworm outbreaks and their qualitative analyses
The ODE model developed in Ludwig et al. (1978) for the dynamics of the budwormforest interaction is a three-variable system based upon the dynamics of the budworm population size and the state of the forest. The state of the forest is determined by the age and size distribution of the trees, their foliage quantity and their physiological condition. As the periodic budworm outbreaks synchronize the development of the trees, the age distribution in the model from Ludwig et al. (1978) is replaced by a single variable S, which represents the average size of the trees. Qualitatively, S can be thought of as the total surface area of the branches in a stand. Similarly, the condition of the foliage and health of the trees is summarized in a single variable E, which represents the 'energy reserve' of the forest. The size of the budworm population is denoted by B. The three-variable ODE model is given in Equations (1) Notice that all of the equations utilize logistic-type growth. The parameters r B , r S and r E represent the intrinsic growth rates for the dynamic variables, K S , and K E are the carrying capacities for S and E. The carrying capacity for B depends upon the surface area of the trees: K B = K S where, in turn, K depends upon the condition of the foliage, so The parameter T E is the value for E where K = K /2. The negative terms on the right-hand side of Equation (1) reflect the reduction in the budworm population due to predation (Crawford & Jennings, 1989). The use of a Hill function appropriately addresses the fact that the rate of predation is bounded since consumption of budworms by individual avian predators is limited by saturation, and the number of birds is limited by factors such as territorial behaviour. In a similar way, the negative term on the right-hand side of Equation (3) reflects a reduction in the rate of growth of the foliage due to damage caused by the budworm population, where the rate of damage is proportional to the number of budworms per unit surface area B/S. Finally, the coefficient of proportionality is modelled by because the stress on the trees is related to the amount of foliage consumed by the insects. Further details can be found in Ludwig et al. (1978). The budworm can increase its density several hundred folds in a few years. Therefore, in a continuous representation of this process, a characteristic time interval for the budworm is of the order of months. The trees cannot put forth foliage at a comparable rate: a characteristic time interval for trees to completely replace their foliage is on the order of 7-10 years. Moreover, the life span of the trees themselves is between 100 and 150 years, in the absence of budworm, so that their generation time is measured in decades. Thus the budworm is a fast variable, while the forest variables E and S are slow variables.
The qualitative analysis performed in Ludwig et al. (1978) follows the usual path: 1) it first examines the dynamics of the fast variable, assuming constant values for the slow variables and 2) the behaviour of the slow variables is analysed while the fast variable is held at its corresponding equilibrium. It is shown in Ludwig et al. (1978), that the Hillfunction terms E 2 /(T E 2 + E 2 ) and (T E 2 + E 2 )/E 2 in Equations (1) and (3) only impact the qualitative behaviour of the model for values of E near zero. As this term is bounded above and below (away from zero) for larger values of E, in the interest of simplifying the qualitative analyses, it could be replaced by a constant. Without loss of generality, we could then assume that for values of E away from zero the terms E 2 /(T E 2 +E 2 ) and (T E 2 +E 2 )/E 2 are bounded above and below by positive constants. In the considerations below, we will assume that those constants have been absorbed by K B and P in Equations (4) and (5). We refer the reader to the original paper for the details, but we will present some highlights from the analyses next, as we will refer to them in Section 3 where our Boolean models are compared qualitatively to those from the ODE model from Equations (1)-(3).
Analysis of the budworm population density. Assuming that S and E are fixed with a value of E = 0 (slow variables), Equation (1) for B (fast variable) becomes where α = α S. The term g(B) = β B 2 α 2 + B 2 models the predation, with α and β representing, respectively, the scale of budworm densities at which saturation begins to take place, and the saturation level of predation.
In addition to the trivial equilibrium B = 0, the equilibria of Equation (6) satisfy the equation Introducing the scaled budworm density μ = B/α and denoting R = (αr B )/β and Q = K B /α, after some simple calculations, one can obtain the following equation for the non-trivial equilibria for B: Equation (7) provides an easy explanation for the bistable behaviour of the system. The left-hand side represents a line that gives the per-capita growth rate of the scaled budworm population and the right-hand side is the death rate due to predation (for the same scaled variable μ). The two sides of the equation are plotted in Figure 1. Depending on the values of R and Q, Equation (7) can have at least one and at most three equilibria. Notice that even though in this case R and Q are parameters for Equation (7), they are in reality functions of the slow variables E and S. Since R is proportional to α, where α = α S, R is dependent on S (and independent from E). On the other hand, the expressions for Q show that it depends on E but not on S. So, we should interpret Figure 1 to mean that for any fixed values of S and E, the value of B will converge rapidly to one of the stable equilibria (the values of which will be B − = αμ − and B + = αμ + ).
The dynamics of the system can be visualized in Figure 1 by imagining that initially B = B − and R is low. If R is then slowly increased (sliding upward) while keeping Q fixed, the values of B − and B c = αμ c will get closer together quite rapidly, compared to a relatively small increase of B + . At a value for R where B − and B c merge, the lower stable equilibrium will no longer exist, and subsequent increases in R will send the budworm population to the equilibrium state B + . If we now start decreasing R slowly (sliding its value down), the value of B + will be decreasing slowly, even beyond the time where the line first touches the curve Figure 1. The death rate due to predation (curve a) and the growth rate (line b) of the scaled budworm density μ. The two stable equilibria occur at μ − and μ + ; the unstable equilibrium is at μ c . Figure reproduced from Ludwig et al. (1978).
μ/1 + μ 2 (corresponding to a state where B − and B c coincide) and then, as R decreases further, intersects the curve in three points, recreating the distinct equilibria B − and B c . Only when R assumes even lower values and B c and B + coincide will the upper equilibrium be lost and the system be driven to equilibrium B − . The argument can be repeated to show that changing Q produces the same effect. Thus, the budworm population is bistable and exhibits the hysteresis property, as described in the Introduction: when the intermediate, unstable equilibrium, B c , coalesces with either the upper or the lower equilibrium, the population may either jump from a low value to a high value or vice versa.
As pointed out in Ludwig et al. (1978), the equilibria B − and B + correspond to budworm limitation by predators and food, respectively. As forest conditions improve, budworm growth exceeds the control by predators and an outbreak occurs. On the other hand, if the forest is destroyed far enough, the predators can again regain control. The conditions at which the jumps may occur depend on the relative relation between R and Q. The critical curves where the upper/lower stable equilibrium coalesces with the unstable equilibrium B c are calculated in Ludwig et al. (1978) and presented in Figure 2.
Analysing the slow variables E and S with the fast variable B held at its equilibrium value. As mentioned earlier, the equations for these variables (Equations (2) and (3)) for values of E away from zero reduce to The multiple K E /E in Equation (8) indicates that the surface area may not be able to increase its carrying capacity if the forest foliage is under stress. However, in endemic times, the value of E will be close to its carrying capacity, and K E /E will be close to 1. The  second term on the right-hand side of Equation (9) denotes the decrease in the foliage destroyed by the budworms. The loss is proportional to B/S -the number of budworms per branch. When E = 0, the coefficient of proportionality P can be considered a constant. Setting dS dt = 0, the null-clines for S, from Equation (8), are For E, after setting dE dt = 0, the null-cline is .
The null-cline for S is a straight line through the origin in the SE-plane, which may or may not intersect the curve from Equation (11). For small values of B, the two null-clines intersect, as shown in Figure 3. The intersection point D is a stable equilibrium and C is unstable, with only a pair of trajectories leading into C (the heavy arrows into C in Figure  3) (Ludwig et al., 1978). If E and S start off to the right of the heavy arrows, they will converge (as t → ∞) to D. If they start off to the left of the heavy arrows, their values decrease.
For large values of B however, the two null-clines will not intersect, as shown in Figure  4. The stable equilibrium is lost and, as t → ∞, the populations S and E are decreasing. For Equations (8) and (9), the decrease in E could lead to E taking negative values, which is unrealistic. The presence of the term E 2 /(T E 2 + E 2 ) in Equations (1) and (3) fixes this problem (see Ludwig et al., 1978 for more details).

Preliminaries
In a Boolean model, all model variables are discretized to take values 1 or 0 to denote high or low levels, respectively. A certain threshold of discretization is chosen for each variable (see, e.g. Dimitrova et al., 2010), and values below the threshold are considered 0, while those above the threshold level. This makes Boolean models specifically suited to model systems with variables that exhibit switch-like behaviour.
As in the differential equation models, in the Boolean model there is one equation per variable. A Boolean model of n-variables x 1 , x 2 , . . . , x n consists of n transition Boolean functions (also called update functions or update rules) f x 1 , f x 2 , . . . , f x n describing the dynamical evolution of the model variables, where f x j (x 1 , x 2 , . . . , x n ) : {0, 1} n → {0, 1}, j = 1, 2, . . . , n, are constructed using the operations AND (∧), OR (∨) and NOT (denoted by a horizontal bar over the respective variable/expression). We should note that in Boolean models elimination is implicit -if maintaining a high level for a variable does not follow from its update rule, it is assumed that, in a time step, the value of the that variable will go below the established threshold value and thus change to zero. Equation (12) presents a simple example of a three-variable model.
Boolean models are time-discrete and the dynamical evolution of the model is determined by iterating the transitions defined by the update functions, starting from an initial condition . Substituting x 0 in the right-hand side of Equation (12) produces the state x 1 of the systems at time t = 1, and so on. Various update schedules are possible (see, e.g. Albert & Robeva, 2015) but, in this paper (except for parts of Section 3.4), we use synchronous update: substituting the values of the variables x t = (x t 1 , x t 2 , . . . , x t n ) at time t produces the new state of the system x t+1 = (x t+1 1 x t+1 2 , . . . , x t+1 n ) at time t + 1. For the discussion that follows, we will focus primarily on the steady states (fixed points) of the model for which the update regime doesn't matter. Cycles, however, may not be preserved when switching from one update schedule to another (see Sections 3.4 and 4 below).
For models with a small number of variables, the transitions determined by the update rules for every possible initial state can be visualized by the (directed) state-space transition graph for the model. The set of vertices of this graph is the set of all possible states x = (x 1 , x 2 , . . . , x n ) ∈ {0, 1} n of the system. An edge between two vertices x = (x 1 , x 2 , . . . , x n ) and y = (y 1 , y 2 , . . . , y n ) in the state-space transition graph indicates that the system can transition from x to y in one step; that is, f x j (x 1 , x 2 , . . . , x n ) = y j , j = 1, 2, . . . , n.
If the update rules for some or all variables do not apply with deterministic certainty but, instead, the transitions occur with specified probabilities, the model is stochastic. For deterministic models, each state in the state-space graph transitions to a unique new state, while in the stochastic case, a state can transition to one of several states with certain probabilities .

A Boolean model without bistability
Using Equations (1)-(3), one can draw the dependency graph for the ODE model presented in Figure 5. Notice that the link from B to E is an inhibition; the rest of the links denote activation. The following model is a direct translation from the dependency graph into Boolean transition rules.
The justification for the Boolean rules is as follows:  Note that the Hill coefficient with a value of 2 in the switching function for bird predation (the negative term in the right-hand side of Equation (6)) is critical in the ODE model, as it generates the three equilibrium states for B, leading to bistability. In the Boolean model, we essentially have a Boolean switch for the elimination, which could be thought of as a Hill function with a very high Hill coefficient.
The state-space transition graph of the Boolean model from Equations (13)-(15) is given in Figure 6. There are two fixed points, for both of which the budworm population settles at B = 0. This could happen due to a lack of food (the fixed point (0, 0, 0)) or because the population of budworms is at a very low level due to predation (the fixed point (0, 1, 1)). Notice that there are two different scenarios under which the system settles at the state (0, 0, 0): 1) when the energy reserve (E = 0) is low, this drives S to zero, even when the insect population is low (B = 0), and 2) when the budworm population is high (B = 1), it eventually consumes the leaves and destroys the forest, before driving itself to extinction (B = 0). On the other hand, if the budworm population is low (B = 0) and the state of the forest is healthy (S = 1 and E = 1) the system is in a steady-state (the fixed point (0, 1, 1) in Figure 6).
Let us consider now the fast variable B as a parameter and examine the behaviour of S and E from Equations (14) and (15) for a fixed value of B, as was done in Ludwig et al. (1978) for the ODE model and discussed above. The results are as follows: When B = 1 is fixed: The transition graph obtained from Figure 6 considering the states where B = 1 is given in Figure 7(a). This qualitative behaviour corresponds to the situation depicted in Figure 4: when B is large, the system settles in the state (S, E) = (0, 0). When B = 0 is fixed: The transition graph obtained from Figure 6 considering the states where B = 0 is given in Figure 7(b). This qualitative behaviour corresponds to the scenario depicted in Figure 3: when B is small, the ODE model features a non-trivial stable equilibrium (the stable equilibrium point D in Figure 3),  corresponding to values (S, E) away from zero. The Boolean model also captures this behaviour, which is given by the existence of the fixed point (1,1) when B = 0.

A Boolean model with bistability
As seen in Section 2, the information from Figures 1 and 2 shows that, based on the value of the variables R and Q, the ODE model can exhibit bistability. Since, as we mentioned in Section 2, R is assumed to be proportional to S, the role of R in the Boolean model will be played by S. Similarly, since Q is proportional to K and K is an increasing function of E, the role of Q in the Boolean model will be played by E. Thus, to show bi-stability in the Boolean implementation of the original ODE model from Ludwig et al. (1978), we will need to recreate the conclusion depicted in Figure 2 with S and E replacing R and Q, respectively. More specifically, we'll need to show that when S is medium and E is high, there are two stable equilibria for B: B = 0 and B = 1.
To do this, we will follow a technique that was introduced in both Veliz-Cuba and Stigler (2011) and Hinkelmann and Laubenbacher (2011) to model low, medium and high values for S. This is done by introducing a new Boolean variable, S h , to indicate high values for S and S =1 to stand for at least medium values for S. With this, we'll have the following: S = 0 S h = 0 low value for S S = 1 S h = 0 medium value for S S = 1 S h = 1 high value for S Table 1. Expected outcome (L = low value for a variable, H = high value for a variable) from a Boolean approximation to match the region of bistability depicted in Figure 2.  Table 1. With this information in mind, we propose a modification of the Boolean model from Equations (13)-(15) given by Equations (16)-(19). Notice that just as in the analysis of the ODE model, we are keeping E and S h fixed, but modelling the dependence of S on the fixed values of E and S h .
The justification for the transition rules is as follows: Transition Rule for B: The term B ∧ S ∧ E indicates reproduction. A population larger than zero is needed for the population of budworms to reproduce and stay high. The term S h ∧ E indicates that even if B is low (a very small number of insect organisms), it will reproduce and grow to above threshold levels when there are plentiful resources (high surface area of branches and high foliage). Transition Rule for S: The term S ∧ E indicates reproduction -E needs to be high for S to reproduce and stay high. Or, S h = 1 implies that S = 1 at the next step.
The state-space transition graph resulting from the model in Equations (16)-(19) is presented in Figure 8. Notice that (1) when E is low (E = 0), the system goes to a fixed point with B = 0 regardless of the values of S and S h . (2) when E is high (E = 1) and S is low (S = 0 and S h = 0), the system settles at a fixed point with B = 0. (3) when E is high (E = 1) and S is high (S = 1 and S h = 1), the system settles at a point with B = 1; that is, the budworm population is high. (4) when E is high (E = 1) and S is medium (S = 1 and S h = 0), the system has two fixed points: one with B = 0 and another one with B = 1.
These outcomes from the model are just as expected, based on the values in Table 1. In addition, if we consider the system settled at the fixed point (B, S, S h , E) = (1, 1, 1, 1), which has high value for S, since (S, S h ) = (1, 1), and then intervened to lower that to medium ((S, S h ) = (1, 0)), the new steady state would be (B, S, S h , E) = (0, 1, 0, 1). On the other hand, if we consider the system settled in (B, S, S h , E) = (0, 0, 0, 1), which has low value for S, since (S, S h ) = (0, 0), and then intervened to increase that to medium ((S, S h ) = (1, 0)), the new steady state for the system would be (B, S, S h , E) = (1, 1, 0, 1).  Thus, which of the two steady states the system will settle in for intermediate values of S depends on its history (hysteresis), which establishes its bistability.

A Boolean model capturing periodic outbreaks
Assume now that S h remains fixed, as in Equation (17), but we choose to also incorporate the dynamics of E in the model. This is well-justified since, as described in Section 1, among the two slow variables, the condition of the foliage changes on a much faster scale than the surface area of the branches. We replace Equation (19) with Equation (20) below, based on the following justification: Transition rule for E: For the foliage to be high, the surface area of the forest needs to be at least medium, and there should be no budworms. The state-space graph of the Boolean model with transition functions given by Equations (16)- (18) and (20) is presented in Figure 9.
We note that the cycle of length 4 that appears in the state-space transition graph for this model is in complete agreement with the behaviour implied by the ODE model and depicted in Figures 2, 3 and 4. For low values of the budworm population (at B = B − ), when S h = 1 and E = 0 (that is, (S, E) is below the upper critical curve in Figure 2 but above the level 1/2 dashed line), the value of E will grow, and (S, E) will settle at the stable equilibrium D in Figure 3. Under those conditions, the size of the budworm population is limited mainly by predation. When, however, the value of E increases to a point where (S, E) crosses over the upper critical curve in Figure 2, the budworm population begins a rapid growth towards the stable equilibrium B + . As this unfolds, the null-clines for S and E begin to separate, eventually resulting in a situation similar to that in Figure 4. When the two null-clines are separated, E begins to decline toward zero. If the lower level critical curve in Figure 2 is crossed, then the budworm population B will decline rapidly. This represents a complete cycle of the growth/decline interaction of the budworm-forest system.
Since the fixed points of the system do not depend on the specific update mode selected for the model, it has been possible for us to work until now with the simplest mathematical assumption of synchronous update. However, it is well documented that in Boolean networks some periodic attractors, like the cycle in Figure 9, can be artifacts of the updating schedule (Klemm & Bornholdt, 2005). In other words, the cycle in question may be due solely to the synchronous update of the Boolean transition rules, and may not appear under other, biologically more realistic, assumptions. Thus, care needs to be taken to examine the nature of the observed cycle, evaluate its robustness, and demonstrate that it would persist under alternative update schemes. To do that, we'll employ the stochastic framework SDDS developed in Murrugarra et al. (2012). SDDS models stochasticity at the functional level by introducing two update probabilities p ↑ , p ↓ for each variable. The parameter p ↑ gives the probability of using the update rule when it would produce a positive change, that is, if the update makes the variable increase its value from 0 to 1. Likewise, the parameter p ↓ gives the probability of using the update function when it makes the variable change its value from 1 to 0. Thus a SDDS in n variables is a collection of triplets where f i is the Boolean rule for the variable x i in the deterministic system and p ↑ i , p ↓ i ∈ [0, 1] for all i = 1, . . . , n. SDDS can be represented as a Markov chain by specifying its transition matrix in the following way. For each variable x i , i = 1, . . . , n, the probability of changing its value is given by and the probability of maintaining its current value is given by Then the transition matrix is given by The dynamics of the stochastic system depends on the transition probabilities a ij , which depend on the propensity values and the update functions. For the system given by Equations (16)-(19), the SDDS format is specified as follows: In order to distinguish between fast and slow variables, and in agreement with the methodology developed in Murrugarra et al. (2012), we have assigned propensity values of p ↑ B = 0.8 and p ↓ B = 0.8 for the faster variable B and propensities of 0.2 (for both p ↑ and p ↓ ) for the slow variables S, S h , and E.
The state-space transition graph under SDDS is given in Figure 10. We observe that the cycle in Figure 9 persists under the stochastic update and thus this attractor is robust to changes in the update schedule; in particular, it will be present under asynchronous update as well. As expected, the two fixed points (0, 0, 0, 0) and (0, 1, 0, 1) are also preserved, but now (0, 1, 0, 1) can be reached from several other states.
To obtain a time-series version of Figure 10, stochastic simulations were performed using SDDS with the parameters given in Equations (22)-(25). When S h = 0, Figure 11 shows that the budworm population will eventually diminish due the scarcity of food, even if one starts with a high population of budworms; the system settles at one of the fixed points depicted in the connected component on the left-hand side in Figure 10. All simulations in Figure 11 were initialized at (1,1,0,1), where is B is high. When on the other hand, we use (0,0,1,0) as the initial state for the simulations, the time series in Figure 12 shows periodic outbreaks that correspond to the cycle in Figure 10.

Discussion
In the last decade, the use of Boolean network models has gained rapid popularity as a useful tool for modelling in molecular biology. In particular, Boolean approximations of ODE models have been used successfully to model gene regulation (Hinkelmann & Laubenbacher, 2011;Robeva & Yildirim, 2013). As this approach has only been applied at the molecular level, we wondered if the methodology could be used at higher levels of biological organization. In particular, we wondered if Boolean models could be successfully utilized to examine the qualitative behaviour of systems that feature populations with levels that are either very low or very high. The budworm outbreaks in spruce and balsam fir forests present one such example. The Boolean models developed here can be considered Boolean approximations of a classical ODE model of the budworm-forest system from (Ludwig et al., 1978).
Our results show that the qualitative dynamics of the budworm-forest system can be recreated completely from Boolean models, employing analyses that parallel those used to analyse the ODE model. In particular, we have shown that the dynamic behaviour of the slow and fast variables for both models is identical, as is the overall dynamic behaviour of the models. Our Boolean models were also shown to correctly capture the system's bistability, as well as its behaviour of periodic insect outbreaks.
One shortcoming of the synchronous Boolean approach, especially when applied to systems evolving on a large timescale, is that it is unrealistic to assume that all transitions will occur with deterministic certainty. Some inherent randomness will always be present, and one way of incorporating this into the model is by introducing propensities for the transitions, as was done in Murrugarra et al. (2012). Using this method, we have performed a verification that the cycle in Figure 9 for the model given by Equations (16)-(18) and (20) is robust to changes in the updating schedule (Murrugarra et al., 2012).
In closing, this work shows that Boolean models can be used successfully at the population level for systems that exhibit switch-like behaviour. In comparison with the ODE models, the Boolean models are more intuitive and can be generalized more easily for larger systems based only on knowledge about the topology of the wiring diagram. We should underscore that although the Boolean approach provides a fairy general framework for modelling dynamical systems, our goal here was to replicate the findings of the specific classical model from (Ludwig et al., 1978). Many more interesting questions could be considered along those lines. For instance, it will be interesting to examine Boolean models that account for the effects of higher trophic levels on the budworm-forest system not included in Ludwig et al. (1978), as well as investigate the value of Boolean models in addressing the broader question of large shifts in a wide range of ecosystems that could be attributed to alternative stable states (Scheffer & Carpenter, 2003;Folke et al., 2004;Biggs et al., 2009). We also feel that the use of the SDDS framework from Murrugarra et al., (2012) that has already been applied successfully to model the bistable dynamics of the phage-λ infection of bacteria and the oscillatory patterns in the p53-mdm2 system hold further potential to make Boolean models at the population level more realistic. However, additional work accounting for the specific biology in a wide spectrum of other systems at different levels of biological organization would be needed in order to adequately assess the broader validity of this approach.