Kinetic Theory Modeling and Efficient Numerical Simulation of Gene Regulatory Networks Based on Qualitative Descriptions

Abstract

In this work, we begin by considering the qualitative modeling of biological regulatory systems using process hitting, from which we define its probabilistic counterpart by considering the chemical master equation within a kinetic theory framework. The last equation is efficiently solved by considering a separated representation within the proper generalized decomposition framework that allows circumventing the so-called curse of dimensionality. Finally, model parameters can be added as extra-coordinates in order to obtain a parametric solution of the model.

1. Introduction

Using mathematical modeling to address large-scale problems in the world of biological regulatory networks has become increasingly necessary given the quantity of data made available by improved technology. In the most general sense, modeling approaches can be thought of as being either quantitative or qualitative. Quantitative methods, such as ordinary differential equations or the chemical master equation, are widespread in the literature [1–11]; when the model is well developed, the detail therein can be incredibly informative. However, these methods are not well suited for all applications. Quantitative models require in-depth knowledge of the reaction kinetics and generally fail as the problem size grows. The alternative approach, qualitative models, does not possess the same amount of detail, but captures the essential dynamics of the system. In addition, qualitative models have a variety of analysis tools that can be applied regardless of the problem size. Gene regulation, as a sub-genre of biological regulatory networks, is characterized by large numbers of interconnected species whose influences depend on passing some threshold, thus largely sigmoidal behaviors. The application of qualitative methods to these systems can be highly advantageous to the modeler.

In this work, we begin by considering the qualitative framework of process hitting, revisited briefly in Section 2.1. A highly flexible model, process hitting captures the most important dynamics of the system with a relatively simple syntax. The structure of this syntax lends itself to powerful static analysis tools, which can be used to answer some of the most important questions about the model, such as steady states or reachability, without constructing the state space. Realistic models in gene regulation are immense and highly interconnected: even when considering a Boolean space, the very enumeration of the possible states of the resulting system creates a combinatorial explosion. This is a frequent obstacle in the field of computer science and has been dubbed the curse of dimensionality. However, there are some questions for which one must access the underlying probability distribution associated with the Markov transitions of the qualitative model as described in Section 2.2. In addition, gaining access to the probability distribution allows for a qualitative and intuitive analysis of the system as a whole. The most pervasive methods have historically been simulation-based, although there are some instances in which this becomes computationally infeasible. Here, we propose a method to solve the system by treating the Markov equations of a process hitting model with numerical techniques. A reduced-basis method, the proper generalized decomposition (PGD), can be used to overcome the curse of dimensionality and provide fast, computationally inexpensive solutions to an otherwise intractable problem, as discussed in Section 2.3. The fact of using numerical approaches usually employed in quantitative analysis to perform qualitative analysis (addressed in general by using discrete techniques, like process hitting) constitutes the main originality and novelty of the present work. The proposed approach allows calculating a sort of qualitative probability distribution that cannot be easily obtained by using standard discrete strategies. In addition, PGD has certain qualities particularly favorable for applications to gene regulatory networks. Unknown parameters can easily be incorporated into the model at the cost of another dimension, as demonstrated in Section 3.

2. Methodologies

2.1. Qualitative Modeling: Process Hitting

Process hitting is a powerful, simple tool for the analysis of large regulatory networks. Historically related to the discrete models of Stuart Kauffman [12] and René Thomas [13], process hitting attempts to address problems of scalability in classical modeling methods, while maintaining the highest degree of expressiveness possible. Formally a subclass of asynchronous automata, it relies on large degrees of abstraction to describe the system as a whole. All interacting species, whether they be enzymes, genes or proteins, are abstracted as sorts. These sorts are then subdivided into processes, which could represent concentration levels, spatial configuration or any other form that has a distinct qualitative impact on the system.

Processes interact with one another via actions, in which processes hit one another to create a bounce to some new level of the same sort at a given rate. For gene regulatory networks, processes are often abstractions of relevant concentration ranges, discretized domains of real numbers and actions representing activation and inhibition reactions. Figure 1 illustrates how to define sorts, processes and actions from a biological understanding of an interaction. Process hitting relies on the initial construction of the most permissive dynamics, otherwise called generalized dynamics, in which no restrictions are placed on the potential behaviors. An example of this can be seen in Figure 1.

The general dynamics may then be successively enriched by the addition of cooperative sorts in order to best capture some known biological behaviors or eliminate undesirable behaviors. Cooperative sorts represent not species, but rather, the combined effects when multiple regulators interact cooperatively on a single target. These sorts are the combined space of the original species; thus, they must be updated, such that the current state of the cooperative sort is compatible with the current state of each of its components. A visual explanation of the construction of a cooperative sort and its refinement of a process hitting model can be found in Figure 2.

Although this is a very simplistic representation of the inner kinetics of a biological process, process hitting semantics allow us to easily model interactions with only partial knowledge of the logical functions encoded therein and pave the way for powerful static analysis techniques in order to study fixed points, reachability and cut sets, which determine the minimum criteria for reachability, in spite of the present combinatorial explosion [15,16]. Examples of process hitting at work can be found in Section 3, where we use static analysis to compare the fitness of the generalized dynamics model with that of the refined model. Furthermore, these tools are freely available online in a software called PINT. We will not attempt to expound completely on the details of process hitting here, but rather, point those interested towards [16] for a formal and thorough introduction to the modeling framework. As we progress to a biological application in Section 3, greater clarity will be given to the concepts described above, including the relevance of cooperative sorts and the power of static analysis. For a more detailed overview of process hitting (PH), the reader can refer to Appendix A.

2.2. Treating Qualitative Systems with Numerical Techniques

In order to address process hitting’s global results, that is the full and complete description of the systems behavior given an initial condition, we must consider the framework in a stochastic context. Process hitting actions move the system from one state to another.

We consider n different chemical species (sorts, including cooperative sorts) ${\mathcal{S}}_i$, i = 1,…, n, each one having K_i + 1 possible levels (processes) s_i ∈ (0, 1, …, K_i), i = 1, …, n, and a set of m reactions (hits) R_j, j = 1, …, m, with propensity a_j encoding the reaction rate:

$${\chi }_{j,1}^{-}{s}_1^{-}{\mathcal{S}}_1+\dots +{\chi }_{j,n}^{-}{s}_n^{-}{\mathcal{S}}_n\stackrel{a_j}{\to }{\chi }_{j,1}^{+}{s}_1^{+}{\mathcal{S}}_1+\dots +{\chi }_{j,n}^{+}{s}_n^{+}{\mathcal{S}}_n$$

(1)

with χ_j,i controlling the appearance of sort S_i in reaction j. For that purpose, χ is a Boolean variable, i.e., χ_j,i = (0, 1).

The system state is defined from z = (s₁, …, s_n). Thus, reaction j transforms the state $\widehat{\mathbf{z}}=(s_1^{-},\dots ,s_n^{-})$ into z, with $\mathbf{z}=\widehat{\mathbf{z}}+{\mathbf{v}}_j$, ${\mathbf{v}}_j=(s_1^{+}-s_1^{-},\dots ,s_n^{+}-s_n^{-})$.

As a memoryless random walk, each action corresponds to a Markov equation, the so-called chemical master equation (CME), which tracks the net change in the probability of existing at a certain state and time Ψ(z, t):

$$\frac{d\mathrm{\Psi }(\mathbf{z},t)}{dt}={\displaystyle \sum _{j=1}^m{a}_j(\mathbf{z}-{\mathbf{v}}_j,t)\mathrm{\Psi }(\mathbf{z}-{\mathbf{v}}_j,t)-a_j(\mathbf{z},t)\mathrm{\Psi }(\mathbf{z},t)}$$

(2)

with the propensity a_j depending on the system state and time.

The result is a system of linear, time-dependent, differential equations, defined given an initial condition. Some of the most famous and broadly-used techniques for addressing problems such as these have been simulation based. Simulation can become computationally expensive with respect to computing time and available memory when addressing highly multidimensional models, as in the case of Equation (2), when the number of species increases. An alternative approach is the direct application of a numerical method to the Markov equations. Here, we propose proper generalized decomposition (PGD) as an effective and well-suited technique for gene regulatory networks.

PGD [10,17–20] is a multi-linear numerical solver that assumes that the target, in this case the probability distribution Ψ(z, t), can be written as a sum of a product of separable functions of the interacting species, F_i(s_i), i = 1, … n, and time, F_t(t):

$$\mathrm{\Psi }(\mathbf{z},t)=\mathrm{\Psi }(s_1,\dots ,s_n,t)\approx {\displaystyle \sum _{k=1}^Q{\psi }^k(\mathbf{z})\cdot F_t^k(t)={\displaystyle \sum _{k=1}^Q{F}_1^k(s_1)\cdot F_2^k(s_2)\cdot \dots \cdot F_n^k(s_n)}}\cdot F_t^k(t)$$

(3)

In essence, PGD is very close to tensor approximation-based strategies that were successfully considered for solving the CME in [21].

PGD is performed iteratively, starting at some arbitrary guess and searching for sets of functions, one functional product at a time, which will minimize the residual of the running sum. These functions are colloquially called modes; however, since the only objective is the reduction of the residual, there is no underlying notion that they represent the greatest source of variance, as is the case with principal component analysis. Although the accuracy increases with the number of modes considered in the finite sum decomposition, we assume that only a limited numberQof functional products are needed to capture the behavior of the system.

The calculation of functions involved in the separated representation is performed on the fly, from the weak form related to the CME (2):

$$\begin{array}{l}{\displaystyle {\int }_{\mathrm{\Omega }\times ℐ}{\mathrm{\Psi }}^{\ast }(\mathbf{z},t)\frac{d\mathrm{\Psi }(\mathbf{z},t)}{dt}d\mathbf{z}dt}\\ ={\displaystyle {\int }_{\mathrm{\Omega }\times ℐ}{\mathrm{\Psi }}^{\ast }(\mathbf{z},t)}\left\{{\displaystyle \sum _{j=1}^m{a}_j(\mathbf{z}-{\mathbf{v}}_j,t)\mathrm{\Psi }(\mathbf{z}-{\upsilon }_j,t)-a_j(\mathbf{z},t)\mathrm{\Psi }(\mathbf{z},t)}\right\}d\mathbf{z}dt\end{array}$$

(4)

where z ∈ Ω = ω₁ × … × ω_n, ω_i = (0, 1, …, K_i) and $t\in \mathrm{I}=(0,\mathcal{T}]$. In fact, the domains ω_i in which the species levels are defined being discrete, the integral in ω_i reduces to a discrete sum.

As just indicated, the algorithm proceeds iteratively, by computing one term of the finite sum at each iteration. If we assume that at iteration p − 1 the solution writes:

$${\mathrm{\Psi }}^{p-1}(\mathbf{z},t)={\displaystyle \sum _{k=1}^{p-1}{\psi }^k(\mathbf{z})\cdot F_t^k(t)={\displaystyle \sum _{k=1}^{p-1}{F}_1^k(s_1)\cdot F_2^k(s_2)\cdot \dots \cdot F_n^k(s_n)\cdot F_t^k(t)}}$$

(5)

at the next iteration p, the solution Ψ^p(z, t) reads:

$$\begin{array}{l}{\mathrm{\Psi }}^p(\mathbf{z},t)={\displaystyle \sum _{k=1}^p{F}_1^k(s_1)\cdot F_2^k(s_2)}\cdot \dots \cdot F_n^k(s_n)\cdot F_t^k(t)\\ ={\displaystyle \sum _{k=1}^{p-1}{F}_1^k(s_1)\cdot F_2^k(s_2)}\cdot \dots \cdot F_n^k(s_n)\cdot F_t^k(t)+F_1^p(s_1)\cdot F_2^p(s_2)\cdot \dots \cdot F_n^p(s_n)\cdot F_t^p(t)\\ ={\displaystyle \sum _{k=1}^{p-1}{\psi }^k(\mathbf{z})}\cdot F_t^k(t)+{\psi }^p(\mathbf{z})\cdot F_t^p(t)\\ ={\mathrm{\Psi }}^{p-1}(\mathbf{z},t)+{\psi }^p(\mathbf{z})\cdot F_t^p(t)\end{array}$$

(6)

where the unknown functions $F_i^p(s_i)$ and $F_t^p(t)$ must be calculated.

The resulting weak form at the present iteration reads:

$$\begin{array}{l}{\displaystyle {\int }_{\mathrm{\Omega }\times ℐ}{\mathrm{\Psi }}^{\ast }(\mathbf{z},t)\left\{F_1^p(s_1)\cdot \dots \cdot F_n^p(s_n)\frac{d{F}_t^p(t)}{dt}\right\}d\mathbf{z}}\\ -{\displaystyle {\int }_{\mathrm{\Omega }\times ℐ}{\mathrm{\Psi }}^{\ast }(\mathbf{z},t)\left\{{\displaystyle \sum _{j=1}^m{a}_j(\mathbf{z}-{\mathbf{v}}_j,t){\psi }^p(\mathbf{z}-{\mathbf{v}}_j)F_t^p(t)-a_j(\mathbf{z},t){\psi }^p(\mathbf{z})F_t^p(t)}\right\}d\mathbf{z}dt}\\ =-{\displaystyle {\int }_{\mathrm{\Omega }\times ℐ}{\mathrm{\Psi }}^{\ast }(\mathbf{z},t)\left\{{\displaystyle \sum _{k=1}^{p-1}{F}_1^k(s_1)\cdot \dots \cdot F_n^k(s_n)\frac{d{F}_t^k(t)}{dt}}\right\}d\mathbf{z}}dt\\ +{\displaystyle {\int }_{\mathrm{\Omega }\times ℐ}{\mathrm{\Psi }}^{\ast }(\mathbf{z},t)\left\{{\displaystyle \sum _{j=1}^m{a}_j(\mathbf{z}-{\mathbf{v}}_j,t){\displaystyle \sum _{k=1}^{p-1}{\psi }^k(\mathbf{z}-{\mathbf{v}}_j)F_t^k(t)-a_j(\mathbf{z},t)}{\displaystyle \sum _{k=1}^{p-1}{\psi }^k(\mathbf{z})F_t^k(t)}}\right\}d\mathbf{z}dt}\end{array}$$

(7)

with:

$$\begin{array}{l}{\psi }^k(\mathbf{z}-{\mathbf{v}}_j)=F_1^k(s_1-{\upsilon }_{j,1})\cdot \dots \cdot F_n^k(s_n-{\upsilon }_{j,n})\hfill \\ {\psi }^k(\mathbf{z})=F_1^k(s_1)\dots F_n^k(s_n)\hfill \end{array}\}k=1,\dots ,p$$

(8)

The simplest test function Ψ^*(z, t) within a Galerkin framework writes:

$$\begin{array}{l}{\mathrm{\Psi }}^{\ast }(\mathbf{z},t)=F_1^{\ast }(s_1)\cdot F_2^p(s_2)\cdot \dots \cdot F_n^p(s_n)\cdot F_t^p(t)+F_1^p(s_1)\cdot F_2^{\ast }(s_2)\cdot \dots \cdot F_n^p(s_n)\cdot F_t^p(t)+\dots \\ +F_1^p(s_1)\cdot F_2^p(s_2)\cdot \dots \cdot F_n^{\ast }(s_n)\cdot F_t^p(t)+F_t^p(s_1)\cdot F_2^p(s_2)\cdot \dots \cdot F_n^p(s_n)\cdot F_t^{\ast }(t)\end{array}$$

(9)

Introducing the test function (9) into the weak form (7) results in a nonlinear problem whose solution requires an appropriate linearization strategy. The simplest strategy consists of an alternated directions fixed point algorithm, which, considering known functions $F_2^p(s_2),\dots ,F_n^p(s_n)$, $F_t^p(t)$ (initialized randomly), calculates $F_1^p(s_1)$. From the just updated $F_1^p(s_1)$ and $F_3^p(s_2),\dots ,F_n^p(s_n)$, $F_t^p(t)$, function $F_2^p(s_2)$ is updated. Then, function $F_3^p(s_3)$ is updated, and so on, until calculating $F_t^p(t)$ from the just updated $F_1^p(s_1),\dots ,F_n^p(s_n)$. Then, the iteration continues until reaching convergence, that is the fixed point that results in functions $F_1^p(s_1),\dots ,F_t^p(t)$. Because the calculation of functions $F_i^p(s_i)$, for i = 1, ⋯, n, follows the same rationale, in Appendix B, we summarize the procedure for calculating those functions $F_i^p(s_i)$, as well as the function depending on time $F_t^p(t)$.

Having calculated Ψ^p(z, t), the enrichment procedure continues for calculating using the same rationale ${\mathrm{\Psi }}^{p+1}(\mathbf{z},t)={\mathrm{\Psi }}^p(\mathbf{z},t)+F_1^{p+1}(s_1)\cdot \dots \cdot F_n^{p+1}(s_n)\cdot F_t^{p+1}(t)$. The enrichment stops as soon as ‖Ψ^Q − Ψ^Q−1‖ < ϵ. Alternative goal-oriented stopping criteria exist and were successfully applied in our former works [22].

For more details on the fundamental and technical details on the PGD, the reader can refer to [23–26] and the references therein.

When considering a network of n species with K possible levels (for the sake of simplicity, we assume for a while that K_i = K, ∀i), solved at P time steps, the resulting complexity scales with Q · (n · K) (the complexity related to the solution of the time problem being negligible) instead of Kⁿ involved in an hypothetical mesh or grid that should be solved P times. The separated representation is sketched in Figure 3, which emphasizes the fact that a 3D problem can be solved as a sequence of one-dimensional problems. If, for a while, we consider K = 3 and 100 terms in the finite sum decomposition, Q = 100, the complexity of the PGD solver will scale as 300n, with n the number of species involved in the network, while standard mesh-based discretizations will scale as 3ⁿ. We can notice that for n ≈ 7, both complexities become equivalent, and for n ≈ 15, the one related to the PGD is three orders of magnitude lower than the one involved in standard mesh-based discretizations.

The PGD approximation converges to the numerical solution of the problem when increasing the number of terms. To understand this fact, it suffices to realize that standard numerical solutions consist of a polynomial up to a certain degree, and a polynomial in many dimensions is no more than a finite sum of terms, each one involving the product of a function of each coordinate. This fact proves the generality of the separated representation involved in the PGD approach. The main difference with respect to standard procedures is that when considering the PGD approach, the different functions involved in these products are assumed unknown and are calculated on the fly.

It is important to note that the probability distribution obtained when applying the PGD constructor converges towards the one related to Equations (1) and (2). If these (Equations (1) and (2)) describe the true probability distribution, then the solution (3) obtained by invoking the PGD constructor will approach the true probability distribution if we consider an adequate number of terms Q in (3).

In what follows, when considering the chemical master equation formalism for qualitative modeling purposes, the possible levels of each species reduce in general to a few states, e.g., (0, 1) for expressing the presence or not of the considered species or (0, 1, 2) for indicating little, moderate or a high number of individuals in the involved species. Within the qualitative modeling framework considered in the present work, a constant rate is associated with each reaction [9]. Each reaction implies given species with given concentration levels (different concentration levels of the same species could define different reactions with their own rates). The considered rates have a significant impact in the systems dynamics. When considering the PGD framework, one can consider these rates as extra-coordinates for calculating all of the dynamics related to each choice of these rates. This possibility is specific to PGD approaches and could be used precisely for better identifying rates from experiments.

In our former works [10], we considered complex functional dependences of kinetic rates, as is usual when addressing quantitative models within the chemical master equation framework; however, and without loss of generality, in the present applications, all of the involved rates are considered constant, as is usual when considering qualitative models. Qualitative descriptions even when addressed with numerical techniques as the one here proposed can describe the main system features in a qualitative way. Analyses requiring a quantitative approach require finer descriptions and determination of the kinetic rates. Each approach (qualitative and quantitative) has its respective domains of applicability, whose comparison is out of the scope of the present work.

3. Application to a Biological Network

It is easier to understand the concepts of process hitting and PGD, as well as to see their individual and combined benefits when seen in action in the context of a realistic application towards a gene regulatory network. Here, we investigate a medium-scale model of the ErbB signaling pathway, which regulates the cell transition from the G1 to S life phase, an important checkpoint that determines whether a cell should divide, delay division or enter a quiescent state. Overexpression of ErbB is associated with many kinds of cancer, and drugs that target it and its receptor are common treatments for breast, lung and colon cancers. The directed graph for this network was taken from [27], where twenty species interact according to Boolean rules. The directed graph can be found in the Appendix for reference. We begin our application by constructing a process hitting model from this Boolean predecessor, taking the most permissive, generalized dynamics, followed by its refinement via the incorporation of cooperative sorts. The impact that this refinement has on both the statistical analysis and application of PGD will be investigated, both in terms of expressiveness and complexity. Finally, the potential of PGD’s capacity to easily incorporate model parameters as extra coordinates will be demonstrated by taking many potential values for the rates of two reactions in the directed graph.

The translation of a Boolean model to the generalized dynamics of process hitting is relatively straight forward, as shown in Figure 2: the absence of an activator effectively serves as an inhibitor and vice versa. The formal relationship between Boolean networks and process hitting can be found in [28]. At this point, we would like to investigate the model to see if it adequately reflects our biological understanding of the system as a whole: are experimentally demonstrated states reachable; are impossible states unreachable; and are there fixed points if steady-state behaviors exist? These questions constitute sanity checks, making sure our model is not essentially flawed from the beginning. The structure of the system (see Appendix C) suggests two species of experimental interest: EGF as an input, having no predecessor, and pRB as output, having no successor. Using these two species, we can easily formulate simple reachability criteria in order to perform sanity checks on our model. We consider a system at rest, in which all components begin in their inactive state. If no changes are made on the input protein, EGF when it is inactive, we expect that the system will remain at rest and that no change is to occur in the output protein. However, if EFG is introduced, the signal should be able to propagate to the output, pRB. In order to be a feasible model, the system must pass these two criteria. Results from static analysis, shown in

Table 1. Results for ErbB models using generalized dynamics and a refinement with cooperative sorts. Here, the two models were tested using three sanity checks related to our biological understanding of the system: the presence of fixed points, the lack of impossible behaviors and the presence of demonstrated behaviors. In order to be considered a functioning model, pRB should remain at rest when the system is universally inactive, including the absence of input protein EGF. However, in the presence of EGF, a signal should be able to propagate through the system, potentially activating pRB. We see that, while the generalized dynamics were able to propagate a signal from EGF to pRB (EGF present), it was not able to prevent sporadic activation of pRB in a system at rest (EGF absent), nor find any fixed points.

Model	Fixed points	EGF absent	EGF present
Gen. Dynam.	0	Fail	Pass
Refined	3	Pass	Pass

, provide good evidence that the generalized dynamics are too permissive and do not accurately capture the biological behaviors, which are essential for a functioning model: not only are we unable to find any fixed points within the system, for which we do expect to find at least one, but the protein pRB may become sporadically activated in a globally inactive system, failing the first sanity check. Therefore, we must refine the model, incorporating the suggested logical gate rules from [27] via cooperative sorts. In doing so, we recapture these vital phenomena, finding three fixed points and passing both sanity checks. These results were obtained in a matter of seconds, using simple commands in freely available software, allowing us to efficiently alter our model before investing time in more computationally expensive analysis.

3.1. Treating Qualitative Systems with Numerical Techniques

The Markov equations of the process hitting actions provide a system of DEs to which we can apply PGD. Each species occupies a dimension of the state space. With two processes to each sort, the final problem is of size 2²⁰, or over one million possible states. The underlying probability distribution is a function of these species and time. Our goal is to approximate this solution by a summation of separable functions:

$$\mathrm{\Psi }(\mathbf{z},t)\approx {\displaystyle \sum _{k=1}^Q{F}_{EGF}^k(s_{EGF})\cdot \dots }\cdot F_{pRB}^k(s_{pRB})\cdot F_t^k(t)$$

(10)

In the case of process hitting containing only the generalized dynamics, this is an appropriate and accurate method. However, once cooperative sorts are incorporated into the qualitative model, the cooperating species can no longer be represented by separable functions. To satisfy the enriched model, we may simply combine those dimensions that participate in cooperative sorts. While this does create vectors that grow exponentially with each added species, it is biologically implausible that more than three or four species would participate in a cooperative influence on a single target. Therefore, we can expect this growth to be cut short long before the dimension of a cooperative sort becomes too large. As we combine the state spaces so that they reflect their cooperative sorts, the error associated with PGD solutions as compared to the solution obtained from simulation techniques decreases.

However, what is to be done when one species participates in multiple cooperative sorts? There are two possibilities: either group species into a macro-species or leave them separated. The former possibility could blow up the state space of the combined species, due to chains of regulations forcing many species to be glued together. Hence, it may result in being unfeasible. The second one does not blow up (the state of individual species remains the same), but leads in general to more terms in the separated representation, due to the need to correct for the correlation effects of the regulation.

The solutions that we obtain from PGD are approximations of the full probability distribution corresponding to the Markov equations created by process hitting. From these probability distributions, we are able to make fast analysis of the global behaviors of the system: rather than being limited to asking the questions answerable using static analysis, a modeler can watch the system evolve through time and make general statements on the qualitative behavior.

3.2. Incorporation of Unknown Parameters

It is often the case, especially in a growing field such as genomics, that elements of a regulatory network are disputed or unknown. Researchers may come to very different conclusions about the parameters that fit a particular system. With simulation techniques, each new set of parameters requires a full repetition of all of the trials, limiting the modeler and leading to ad hoc choices made for the sake of feasibility. However, PGD offers a simple way of incorporating these unknown parameters directly into the model, making it possible to obtain an approximate solution for a range of values all at once [10,29]. The parameter is encoded as one of the separable spaces and is included at the cost of one dimension added to the overall solution space. For our example, perhaps one of the regulating reactions is difficult to study separately from the system as a whole, say interactions involving p27 and p21. Unlike the first half of the directed graph, which is simply an activation cascade, these proteins are involved in both inhibiting and activating relationships, so changes to their rate laws should more greatly influence the final expression of pRB. We would like to incorporate many potential values of the action firing rate r into our model, anywhere between two times faster 2r and two times slower 2/r than the other reactions in the system. Since our representation requires discretization, we consider forty equally-spaced values between r/2 and 2r. Our decomposition of Ψ(z, t) is changed slightly in order to accommodate the parameter r for the range of possible values from Ψ(z, r, t):

$$\mathrm{\Psi }(\mathbf{z},r,t)\approx {\displaystyle \sum _{k=1}^Q{F}_{EGF}^k(s_{EGF})}\cdot \dots \cdot F_{pRB}^k(s_{pRB})\cdot F_r^k(r)\cdot F_t^k(t)$$

(11)

While simulation run time grows linearly with each element, 40-times longer since there are 40 values in the discretization, to obtain a result, we are able to derive a solution in relatively equal time using PGD. In Figure 4, we see three solutions for the protein pRB given different values of the rate parameter r: r/2, 3r/2 and 2r. These comparably fast results allow us to perform general analysis on the network by directly observing how the global behavior changes with parameters.

4. Evaluation and Analysis

Up to now, we have presented a new method of approaching discrete models of gene regulatory networks, uncovering briefly the origins of its individual components, process hitting and PGD, and applying it to a real biological system. As bioinformatics grows and many new methodologies are proposed, however, it is increasingly important to discuss in a straightforward manner how well our method performs, highlighting both its merits and weaknesses.

We began our approach by considering the discrete modeling framework of process hitting. This approach allows for a great extent of abstraction in the development of a regulatory model and comes with a well-developed analysis toolbox, making it an attractive starting framework. However, the application of PGD to Markov equations would be effective for any discrete modeling type that can be described as such. Process hitting does have an advantage in that the species that interact in nonlinear ways and thus must be represented together in the decomposition are well defined as cooperative sorts in the very construction of the model. Once the Markov equations have been provided, the method is relatively straightforward; PGD is a well-founded numerical method with thoroughly documented implementations.

As for the results themselves, there are several points to touch on: ease of analysis, model validation and accuracy. One of the most interesting aspects of this approach is the nature of the results: a full probability distribution as an approximate solution to a set of equations. As demonstrated in Figure 4, the output lends itself to visualization on an individual scale, that is for each protein involved. The behavior of a gene or protein can be described in very plain and qualitative terms, even for elements whose evolution is complicated and never reaches steady-state behaviors. However, the apparent behavior can only be taken at face value: while a protein may appear to oscillate or tend to a certain value, the solution is only valid for the limited time frame in which it is analyzed and is not tied to a mathematical principle governing its evolution. Furthermore, since the solution is approximative, it is possible that species whose static analysis proves total inactivation or activation would be found slightly activated or inactivated in the PGD solution, a potentially important distinction in the world of genomics. In such instances, the static analysis and numerical analysis may be found in conflict with one another, compromising the accuracy of the resulting analysis.

5. Conclusions

In the case of gene regulatory networks, there are many reasons why a modeler might choose the application of qualitative methods, one of which is process hitting. Process hitting offers many advantages for large-scale, which are often the more realistic, systems in the form of static analysis tools. These analysis tools alone, however, cannot provide the complete and intuitive solution of the system as a full probability distribution for each state over time. By translating process hitting actions to Markov equations, we are able to treat a system of PDEs directly. Proper generalized decomposition has proven efficient in solving process hitting models. As opposed to simulation techniques, which have been historically the preferred methodology, PGD can provide full solutions, including multiple unknown parameters, with a single run. Here, we have shown some of the potential of this method, applying a combination of static analysis and numerical tools in order to maximize the expressiveness and understanding of a qualitative model. Only the basic elements of process hitting have been incorporated into the Markov equations considered, that is actions with simple rate laws. Including temporal and varied stochastic features into these equations would further increase its potential.