Piecewise linear and Boolean models of chemical reaction networks

Piecewise linear and Boolean models of chemical reaction networks

Alan Veliz-Cuba111Correspondence Author
Alan Veliz-Cuba (alanavc@math.uh.edu), Ajit Kumar (ajit.kumar@snu.edu.in), Krešimir Josić (josic@math.uh.edu)
Ajit Kumar Krešimir Josić

Models of biochemical networks are frequently high-dimensional and complex. Reduction methods that preserve important dynamical properties are therefore essential in their study. Interactions between the nodes in such networks are frequently modeled using a Hill function, . Reduced ODEs and Boolean networks have been studied extensively when the exponent is large. However, the case of small constant appears in practice, but is not well understood. In this paper we provide a mathematical analysis of this limit, and show that a reduction to a set of piecewise linear ODEs and Boolean networks can be mathematically justified. The piecewise linear systems have closed form solutions that closely track those of the fully nonlinear model. On the other hand, the simpler, Boolean network can be used to study the qualitative behavior of the original system. We justify the reduction using geometric singular perturbation theory and compact convergence, and illustrate the results in networks modeling a genetic switch and a genetic oscillator.

1 Introduction

Accurately describing the behavior of interacting enzymes, proteins, and genes requires spatially extended stochastic models. However, such models are difficult to implement and fit to data. Hence tractable reduced models are frequently used instead. In many popular models of biological networks, a single ODE is used to describe each node, and sigmoidal functions to describe interactions between them. Even such simplified ODEs are typically intractable, as the number of parameters and the potential dynamical complexity make it difficult to analyze the behavior of the system using purely numerical methods. Reduced models that capture the overall dynamics, or allow approximate solutions can be of great help in this situation [Verhulst (2006); Hek (2010)].

Analytical treatments are possible in certain limits. The approaches that have been developed to analyze models of gene interaction networks can be broadly classified into three categories [Polynikis et al. (2009)]: Quasi Steady State Approximations (QSSA), Piecewise Linear Approximations (PLA), and discretization of continuous time ODEs. In particular, in certain limits interactions between network elements become switch–like [Kauffman (1969); Snoussi (1989); Mochizuki (2005); Alon (2006); Mendoza and Xenarios (2006); Davidich and Bornholdt (2008); Wittmann et al. (2009); Franke (2010); Veliz-Cuba et al. (2012)]. For instance, the Hill function, , approaches the Heaviside function, , in the limit of large . In this limit the domain on which the network is modeled is naturally split into subdomains: The threshold, corresponding to the parameter in the Hill function, divides the domain into two subdomains within which the Heaviside function is constant. Within each subdomain a node is either fully expressed, or not expressed at all. When is large, the Hill function, , is approximately constant in each of the subdomains, and boundary layers occur when is close to the threshold,  [Ironi et al. (2011)]. To simplify the system further, we can map values of below the threshold to 0, and the values above the threshold to 1 to obtain a Boolean network (BN); that is, a map

where each function describes how variable qualitatively depends on the other variables [Glass and Kauffman (1973); Snoussi (1989); Thomas (1990); Edwards and Glass. (2000); Edwards et al. (2001)]. Such reduced systems are simpler to analyze, and share the dynamical properties of the original system, if the reduction is done properly.

The reduced models obtained in the limit of a large Hill coefficient, have a long and rich history. Piecewise linear functions of the form proposed in [Glass and Kauffman (1973)] have been shown to be well suited for the modeling of genetic regulatory networks, and can sometimes be justified rigorously [De Jong et al. (2004)]. In particular, singular perturbation theory can be used to obtain reduced equations within each subdomain and the boundary layers, and global approximations within the entire domain [Ironi et al. (2011)]. On the other hand, although BNs have been used to model the dynamics of different biological systems, their relation to more complete models was mostly demonstrated with case studies, heuristically or only for steady states [Glass and Kauffman (1973); Snoussi (1989); Thomas (1990); Albert and Othmer (2003); Mendoza and Xenarios (2006); Davidich and Bornholdt (2008); Abou-Jaoudé et al. (2009, 2010); Wittmann et al. (2009); Franke (2010); Veliz-Cuba et al. (2012)].

Here we again start with the Hill function, , but instead of assuming that is large, we assume that is small. This case has a simple physical interpretation: Consider the Hill function that occurs in the Michaelis-Menten scheme, which models the catalysis of the inactive form of some protein to its active form in the presence of an enzyme. When is small the total enzyme concentration is much smaller than the total protein concentration. Although the subsequent results hold for any fixed , for simplicity we assume .

More precisely, we consider a model biological network where the activity at each of nodes is described by , and evolves according to


where , and the functions , are affine functions.

This type of equations have been used successfully in many models [Goldbeter and Koshland (1981); Goldbeter (1991); Novak et al. (2001); De Jong (2002); Tyson et al. (2003); Ishii et al. (2007); Ciliberto et al. (2007); Davidich and Bornholdt (2008); van Zwieten et al. (2011)]. Here and describe how the other variables affect and can represent activation/phosphorylation/ production and inhibition/dephosphorylation/decay, respectively. The variables can represent species such as protein concentrations, the active form of enzymes, or activation level of genes. A simple example is provided by a protein that can exist in an unmodified form, and a modified form, (e.g. proteases, and Cdc2, Cdc25, Wee1, and Mik1 kinases [Goldbeter (1991); Novak et al. (1998, 2001)]) where the conversion between the two forms is catalyzed by two enzymes, and [Goldbeter and Koshland (1981); Goldbeter (1991); Novak et al. (1998, 2001)] (See Appendix for details). However, note that the models of chemical reactions we consider can be rigorously derived from the Chemical Master Equation only in the case of a single reaction [Kumar and Josić (2011)]. The models of networks of chemical reactions that we take as the starting point of our reduction should therefore be regarded as phenomenological.

It is easy to show that the region is invariant so that Eq. (1) is a system of equations on . Equations involving this special class of Hill functions are generally referred to as Michaelis-Menten type equations, and the Michaelis-Menten constant [Michaelis and Menten (1913); Goldbeter and Koshland (1981); Goldbeter (1991); Novak and Tyson (1993); Novak et al. (2001); Tyson et al. (2003); Ciliberto et al. (2007); Davidich and Bornholdt (2008); Ma et al. (2009)].

The constants are frequently very small in practice [Novak et al. (2001); Davidich and Bornholdt (2008)], which motivates examining Eq. (1) when . In this case, we discuss a two step reduction of the model

We first illustrate this reduction using two standard examples, and then provide a general mathematical justification. We note that the reduction obtained in the first step (see Eq. (14a)) is actually (algebraic) piecewise affine. However, it is customary to refer to the equation and the associated model as piecewise linear [Glass and Kauffman (1973); Snoussi (1989); Thomas (1990); Edwards and Glass. (2000); De Jong (2002)], and we follow this convention.

The main idea behind the piecewise linear (PL) reduction is simple: If then the Hill functions, . However, when and are comparable, this is no longer true. In this boundary layer, we rescale variables by introducing . A similar argument works for the function (see Appendix). We show that using this observation, the domain naturally decomposes into a nested sequence of hypercubes. The dynamics on each hypercube in the sequence is described by a solvable differential-algebraic system of equations. The PL reduction therefore gives an analytically tractable approximate solution to the original system.

In the next step of the reduction we obtain a Boolean Network (BN): The PL approximation is used to divide into chambers. Within nearly all of a chamber the rate of change of each element of the network is constant when . We use these chambers to define a BN. A similar approach was recently used to motivate a Boolean reduction of a model protein interaction network [Davidich and Bornholdt (2008)].

The mathematical justification also follows two steps. We use Geometric Singular Perturbation Theory (GSPT) in Section 4.1 to justify the PL approximation. The justification of the BN reduction is given in Section 4.2. We show that there is a one-to-one correspondence between steady states (equilibrium solutions) of the BN and the full and PL system near the vertices of . Futhermore, we show that this one-to-one correspondence between steady states is actually global (up to a set of small measure in ). BNs have been used to study oscillatory behavior [Li et al. (2004); Abou-Jaoudé et al. (2009)], and we prove in Section 4.2.3 that under some conditions oscillations in a BN correspond to oscillations in the full system.

2 Example problems

We start by demonstrating the main idea of our approach using networks of two and three mutually repressing nodes. These nodes can represent genes that mutually inhibit each other’s production [Gardner et al. (2000); Elowitz and Leibler (2000)]. However, the theory we develop applies whenever the heuristic model given in Eq. (1) is applicable. We accompany these examples with a heuristic explanation of the different steps in the reduction.

Figure 1: (a) Nodes , inhibiting each others activity resulting in a switch. The node which starts out stronger suppresses the activity of the other. (b) Nodes , , and suppress each other in a cyclic fashion. Under certain conditions, this can lead to oscillations.

2.1 A network of two mutually inhibiting elements

We start with the common toggle switch motif, i.e a network of two mutually repressing elements (see Fig. 1a)  [Tyson et al. (2003); Gardner et al. (2000)]. Let represent the normalized levels of activity at the two nodes. Therefore, when the network element is maximally active (expressed). The activity of the two nodes in the system can be modeled by


where is some positive constant. The structure of Eq. (2) implies that the cube is invariant (see Proposition 1).

2.1.1 Piecewise linear approximation

In the limit of small , Eq. (2) can be approximated by a piecewise linear differential equation: If is not too close to zero the expression is approximately unity. More precisely, we fix a small , which will be chosen to depend on . When and is small then . Similarly, when then .

With this convention in mind we break the cube into several subdomains, and define a different reduction of Eq. (2) within each. Let to denote the region where is the set of variables that are close to 0, and is the set of variables close to 1 (See Table 1 and Eq. (3.1)). Also, we omit the curly brackets and commas in (e.g. and ) (see Fig. 2).

Figure 2: Subdomains for the unit square .

We first reduce Eq. (2) on each of the subdomains. The interior of the domain consist of points where neither coordinate is close to 0 nor 1, and is defined by


Eq. (2), restricted to is approximated by the linear differential equation


On the other hand, if one of the coordinates is near the boundary, while the other is in the interior, the approximation is different. For instance, the region


forms a boundary layer where is of the same order as .

The term cannot be approximated by unity. Instead the approximation takes the form


This equation can be simplified further. Since is invariant (for ), must be small inside the boundary layer (see Fig. 4). We therefore use the approximations in Eq. (6a) and in Eq. (6b) to obtain


Note that Eq. (7a) is linear and decoupled from Eq. (7b), while Eq. (7b) is an algebraic system which can be solved to obtain . Within we thus obtain the approximation


We only have the freedom of specifying the initial condition , since is determined by the solution of the algebraic equation (7b). As we explain below, this algebraic equation defines a slow manifold within the subdomain . The reduction assumes that solutions are instantaneously attracted to this manifold.

Table 1 shows how this approach can be extended to all of . There are 9 subdomains of the cube, one corresponding to the interior and four each to the edges and vertices. On the latter eight subdomains, one or both variables are close to either 0 or 1. Following the preceding arguments, variable(s) close to 0 or 1 can be described by an algebraic equation. The resulting algebraic-differential systems are given in the last column of Table 1. Furthermore, by using the approximations for and for , we obtain a simple approximation of the dynamics in each subdomain which is -th order in . For example, in , we obtain the approximation .

Each approximate solution can potentially exit the subdomain within which it is defined if at some time or and the -th coordinate of the vector field is positive or negative, respectively. This can happen when the sign of some entry of the vector field changes; that is, solutions can exit subdomains when they reach a nullcline. Also, solutions can leave the subdomain if they started on the other side of the nullcline to begin with. The global approximate solution of Eq. (2) is obtained by using the exit point from one subdomain as the initial condition for the approximation in the next. In subdomains other than some of the initial conditions will be prescribed by the algebraic part of the reduced system. The global approximation may therefore be discontinuous, as solutions entering a new subdomain are assumed to instantaneously jump to the slow manifold defined by the algebraic part of the reduced system. Fig. 3 shows that when is small, this approach provides a good approximation.

Table 1: List of differential–algebraic systems that approximate Eq. (2) in different parts of the domain. The subdomains are named so that the superscript (subscript) lists the coordinates that are close to (close to 0), with 0 denoting the empty set. For example, denotes that subdomain with and , and the subdomain where is near , but is away from the boundary. The middle column define the subdomain explicitly. The right column gives the differential-algebraic system that approximates Eq. (2) within the given subdomain.
Figure 3: Comparison of the numerical solution of Eq. (2) (dashed black) and the solution of the approximate system as listed in Table 1 (colors) for two different values of . We used in (a); and in (b). Different colors are used for the solution of the reduced system in different subdomains. Solution of the linear approximation started in the subdomain (Initial value: ), and as soon as decreased below , we assumed that the solution entered subdomain . The approximate solution is discontinuous since when , the solution jumped (see inset) to the manifold, described by the algebraic part of the linear differential algebraic system prevalent in the subdomain , Eq. (7b). The solution finally stopped in the subdomain .

2.1.2 Boolean approximation

We now derive a Boolean approximation, , that captures certain qualitative features of Eq. (2). The idea is to project small values of to 0 and large values of to 1, and map the value of the -th variable into 0 and 1 depending on whether is decreasing or increasing, respectively. We will show that the resulting BN can be used directly to detect steady states in the corner subdomains.

Note that for a BN time is discrete; a time step in the Boolean approximation can be interpreted as the time it takes the original system to transition between chambers. Different transitions in the Boolean network may have different duration in the original system; so the time steps in the BN are only used to keep track of the sequence of events, but not their duration.

The reduction described in the previous section gives a linear ODE in the interior region (Eq. (4)), where approaches as . The approximating linear system therefore provides significant information about the behavior of the full, nonlinear system for small.

We first examine the nullclines. In Fig. 4 we can see that as decreases, in the interior of the nullclines of Eq. (2) approach the nullclines of Eq. (4) given by and restricted to . These lines divide the domain into four chambers, which we denote

Figure 4: Behavior of nullclines as decreases. Top: Nullclines of Eq. (2) for (left) and (right). Bottom: Nullcline of Eq. (2) (black curve) and the manifold defined by Eq. (7b) (red) for (left), and (right).

On the other hand, the part of the nullclines inside the boundary subdomains are approximately the slow manifolds defined by equivalents of Eq. (8b). Here the slow manifolds converge to the nullclines as (See Fig. 4).

As a shorthand, we define the “sign” of a vector as the vector composed by the signs of its components, . Note that although the sign of the vector is constant in each chamber, the sign of the vector field of Eq. (2) may differ. For example, in chamber , the sign of the vector field can take all possible values. However, this difference is small when is small, because the regions between the nullclines approach the actual chambers (Fig. 4).

We consider Eq. (2) in each chamber, starting with the first coordinate, . For any solution with initial condition in , the sign of is positive and increases within the chamber. We use this observation to define . The formal definition of this function will be given below – intuitively maps a chamber to 1 if is increasing within the chamber, and to 0 otherwise. Similarly, since initially increases within , we let . Similarly we set , , , and . The -th variable is “discretized,” i.e. mapped to 0 and 1 depending on whether is decreasing or increasing, respectively.

More formally, consider the set , with each element identified with a chamber (e.g., the element represents the chamber ). Then and are defined as Boolean functions from to by setting , , , and , , . These two Boolean functions define a BN, . The functions also define a dynamical system, . However, other update schedules can be used [Aracena et al. (2009)].

The BN reduction can be obtained easily from the sign of at the vertices of , since the sign of the vector field is constant within a chamber. To do so we use the the Heaviside function, , defined by if , if , and . For example, in , both entries increase. We can see this by evaluating for . Using the same argument in each chamber, we obtain the BN


where we used the convention that acts entrywise on each component in the argument.

2.1.3 Steady states of the BN and the PL approximation

While the BN gives information about which variables increase and decrease within a chamber, it is not yet clear how or if the dynamics of the BN in Eq. (9), and the PL approximation in Table 1 are related.

We next show that the steady states of the PL approximation near the vertices can be determined by the steady states of the BN. The reduced equations in the corner subdomains and are purely algebraic. When is small, some of these equations have a solution in , indicating a stable fixed point near the corresponding corner (in this case and ). Others will not have a solution in , indicating that approximate solutions do not enter the corresponding subdomain (here and ). To make the relationship between steady states less dependent on the actual parameters, consider the system

where and . In the previous example , and .

Now, at the corner subdomain we have the approximate equations,

or equivalently,

These equations have a solution in if and only if and , or equivalently, if and only if . A similar analysis leads to the following conditions for the existence of fixed points in each corner subdomain

More compactly, the condition is , where is the corner of interest. Hence, the BN can also be used directly to detect which corner subdomains contain steady states.

The relationship between steady states in the full system at the corner subdomains and the steady states of the BN is straightforward. However, since there are many update schemes for BNs, the relationship between trajectories is more subtle. For example, using synchronous update we obtain the transition which is not compatible with the solutions of Eq. (2) (See Fig. 5). On the other hand, using asynchronous update we obtain the transitions and which are more representative of the solutions of Eq. (2). Thus, we will focus on transitions that are independent of the update scheme, that is, transitions where only one entry changes.

Figure 5: Left: Solutions of Eq. (2) for . When a solution is close to the boundary regions of and , they enter the invariant region as shown in Fig. 3b. Right: Graphical representation of the Boolean transitions (, , , ).

2.2 A network of three mutually inhibiting elements

The same reduction can be applied to systems of arbitrary dimension. As an example consider the repressilator [Tyson et al. (2003); Elowitz and Leibler (2000)] (see Fig. 1a) described by


The cyclic repression of the three elements in this network leads to oscillatory solutions over a large range of values of . The domain of this system, , can be divided into 27 subdomains corresponding to 1 interior, 6 faces, 12 edges, and 8 vertices.

We can again approximate Eq. (2.2) with solvable differential–algebraic equation within each subdomain, to obtain a global approximate solution (See Fig. 6). Note that both the numerically obtained solution to Eq. (2.2) and the solution to the piecewise linear equation exhibit oscillations, and that the approximation is discontinuous. Again, in the limit we obtain a continuous 0-th order approximation.

In this singular limit, solutions can exit a subdomain when they reach a nullcline of the linear system. For example, when is close to 0 and a solution transitions from to , the sign of the second entry of changes from negative to negative; so the second entry of the solution starts increasing (see Fig. 6, panel (e)). Solutions therefore leave the subdomain on which is small and enter the subdomain where . Similarly when is close to 1 solutions transitions from to , and the sign of the first entry of changes from positive to negative. Hence the first entry of the solution starts decreasing (see Fig. 6, panel (f)), and solutions leaves the subdomain where and enter another where .

The BN corresponding to Eq. (2.2), is given by , where is the Heaviside function acting entry wise on the arguments. Eq. (2.2) does not have steady states at the corner subdomains, and neither does the corresponding BN. A subset of states belong to a periodic orbit of the BN:

Note that subsequent states in this orbit differ in a single entry. Thus, the transitions between the states have an unambiguous interpretation in the original system: The BN predicts that if the initial condition is in chamber , then solutions of Eq. (2.2) will go to chamber , then to , and so on. Indeed, solutions of Eq. (2.2), are attracted to a periodic orbit that transitions between the chambers in this order. The remaining two states form a period two orbit under synchronous update, . Here the BN does not give precise information about the dynamics of the original system. We will show that under certain conditions, orbits of the BN where only entry changes at each timestep, correspond to oscillations in the original system.

Figure 6: Comparison of the numerical solution of Eq. (2.2) (dashed black) and the solution of the approximate piecewise linear system (colors) for two different sets of and . For (a)-(c) ; for (e)-(g) . The approximate solution changes color when switching between different subdomains. Panel (d) shows the time series for the solutions of Eq. (2.2) (top) and the PL system (bottom), corresponding to (a)-(c). Panel (h) shows the time series for the solutions of Eq. (2.2) (top) and the PL system (bottom), corresponding to (e)-(g).

3 General reduction of the model system

The approximations described in the previous section can be extended to the more general model given in Eq. (1):

where are some positive constants. Here and are activation/inhibition functions that capture the impact of other variables on the evolution of . The initial conditions are assumed to satisfy for all .

We assume that the activation and inhibition functions are both affine [Novak et al. (2001); De Jong (2002)],


where we use the convention and . The matrix, and the vector capture the connectivity and external input to the network, respectively. In particular, gives the contribution of the variable to the growth rate of variable. If , then appears in the activation function for ; and if then appears in the inhibition function for . The intensity of the external input to the element is , and it contributes to the activation or the inhibition function, depending on whether or , respectively.

Proposition 1.

The cube is invariant for Eq. (1).


It will be enough to show that the vector field at any point on the boundary is not directed outward. Since, and , for any ,

3.1 The PL approximation

To obtain a solvable reduction of Eq. (1) we follow the procedure outlined in Section 2. We first present the results, and provide the mathematical justification in the next section. For notational convenience we let , with . The general case is equivalent. We will use to define the thickness of the boundary layers. We start with the subdivision of the -dimensional cube, .

Let and be two disjoint subsets of , and let


We extend the convention used in Table 1, and in Eqs. (3) and (5) so that when is empty; when is empty; and when , are both empty.

Within each subdomain , Eq. (1) can be approximated by a different linear differential–algebraic system. Following the reduction from Eq. (2) to Eq. (6), for we obtain the linear system

For one of the nonlinear terms remains and we obtain
while for we will have

Eq. (13) is simpler than Eq. (1), but it is not solvable yet. Following the reduction from Eq. (6) to Eq. (7), we now further reduce Eqs.(13b13c). First we use the approximations and in the activation and inhibition functions appearing in Eq. (13). Second, we assume that for and for are in steady state.

Under these assumptions we obtain the reduction of Eq. (1) within any subdomain


Eq. (14) is solvable since Eq. (14a) is decoupled, and Eqs.(14b) and (14c) are solvable for and , respectively, as functions of the solution of Eq. (14a).

Note that in the singular limit we obtain the -th order approximations:

3.2 Boolean approximation

To obtain the Boolean approximation we follow the process described in Section 2. We consider the chambers determined by the complement of the union of the hyperplanes (restricted to ) where . We denote with the set of all chambers ; alternatively, is the set of connected components of . We assume that for all and for all . This guarantees that each corner of belongs to a chamber. The set of parameters excluded by this assumption has measure zero.

Let and be two disjoint subsets of such that and let be the corner that belongs to the corner subdomain . Note that for and for . The chamber that contains the corner in subdomain will be denoted by . We do not name the remaining chambers.

In general, the chambers can be more complex than in the examples of Section 2. Chambers do not have to be hypercubes, different corners may belong to the same chamber, and some chambers may not even contain a corner of , as illustrated in Fig. 7. In the first example, and belong to the same chamber, that is, , and neither containing nor containing are rectangles. Also, has three elements: , , and . In the second example, two chambers do not contain any corner of and are not named. Hence, has four elements: , , and two unnamed chambers that contain no corners.

Figure 7: Chambers for and (left); and (right).

To define the BN, at , we need to find the signs of the components of the vector field on the chamber that contains . Consider . Within the signs of the components of do not change and are equal to the signs of the components of . If the sign of the -th component is negative, we let and if the sign is positive we let . In general, we can write


Hence the value of the BN at a corner is given by the Heaviside function, applied entrywise to . Note that corners that are in the same chamber get mapped to the same point.

Importantly, using Eq. (16) we can compute the BN directly from Eq. (1). For example, for Eq. (2) we have ; and for Eq. (2.2) we have .

Below we show that up to a set of small measure, the BN preserves information about the steady states of the original system. We will also show that under some conditions, “regular” trajectories of a BN correspond to trajectories in the original system.

4 Mathematical justification

We next justify the different approximations made above: In Section 4.1 we use Geometric Singular Perturbation Theory (GSPT) to justify the PL approximation. In Section 4.2 we show that steady state information is preserved by the BN and that, under certain conditions, the BN also provides qualitative information about the global dynamics of the original system.

4.1 Piecewise linear approximation

To obtain the reduced equations at the boundary of , we define the following rescaled variables


Using Eq.  (17) in Eq. (13) we get for

and for ,