Estimating Propensity Parameters Using Google PageRank and Genetic Algorithms

Stochastic Boolean networks, or more generally, stochastic discrete networks, are an important class of computational models for molecular interaction networks. The stochasticity stems from the updating schedule. Standard updating schedules include the synchronous update, where all the nodes are updated at the same time, and the asynchronous update where a random node is updated at each time step. The former produces a deterministic dynamics while the latter a stochastic dynamics. A more general stochastic setting considers propensity parameters for updating each node. Stochastic Discrete Dynamical Systems (SDDS) are a modeling framework that considers two propensity parameters for updating each node and uses one when the update has a positive impact on the variable, that is, when the update causes the variable to increase its value, and uses the other when the update has a negative impact, that is, when the update causes it to decrease its value. This framework offers additional features for simulations but also adds a complexity in parameter estimation of the propensities. This paper presents a method for estimating the propensity parameters for SDDS. The method is based on adding noise to the system using the Google PageRank approach to make the system ergodic and thus guaranteeing the existence of a stationary distribution. Then with the use of a genetic algorithm, the propensity parameters are estimated. Approximation techniques that make the search algorithms efficient are also presented and Matlab/Octave code to test the algorithms are available at http://www.ms.uky.edu/~dmu228/GeneticAlg/Code.html.


Introduction
Mathematical modeling has been widely applied to the study of biological systems with the goal of understanding the important properties of the system and to derive useful predictions about the system.The type of systems of interest ranges from the molecular to ecological systems.At the cellular level, gene regulatory networks (GRN) have been extensively studied to understand the key mechanisms that are relevant for cell function.GRNs represent the intricate relationships among genes, proteins, and other substances that are responsible for the expression levels of mRNA and proteins.The amount of these gene products and their temporal patterns characterize specific cell states or phenotypes [24].
Gene expression is inherently stochastic with randomness in transcription and translation.This stochasticity is usually referred to as noise and it is one of the main drivers of variability [27].Variability has an important role in cellular functions, and it can be beneficial as well as harmful [7,14].Modeling stochasticity is an important problem in systems biology.Different modeling approaches can be found in the literature.Mathematical models can be broadly divided into two classes: continuous, such as systems of differential equations and discrete, such as Boolean networks and their generalizations.This paper will focus on discrete stochastic methods.The Gillespie algorithm [8,9] considers discrete states but continuous time.In this work, we will focus on models where the space as well as the time are discrete variables.For instance, Boolean networks (BNs) are a class of computational models in which genes can only be in one of two states: ON or OFF.BNs and, in general, multistate models, which allow genes to take on more than two states, have been effectively used to model biological systems such as the p53-mdm2 system [1,5,25], the lac operon [37], the yeast cell cycle network [18], the Th regulatory network [21], A. thaliana [3], and many other systems [2,6,10,11,30,39].
Stochasticity in Boolean networks has been studied in different ways.The earliest approach to introduce stochasticity into BNs was the asynchronous update, where a random node is updated at each time step [34].Another approach considers update sequences that can change from step to step [22,29].More sophisticated approaches include Probabilistic Boolean Networks (PBNs) [31] and their variants [17,19].PBNs consider stochasticity at the function level where each node can use multiple functions with a switching probability from step to step.SDDS [25] is a simulation framework similar to PBNs but the key difference is how the transition probabilities are calculated.SDDS considers two propensity parameters for updating each node.These parameters resemble the propensity probabilities in the Gillespie algorithm [8,9].This extension gives a more flexible simulation setup as a generative model but adds the complexity of parameter estimation of the propensity parameters.This paper provides a method for computing the propensity parameters for SDDS.
For completeness, in the following subsection, we will define the stochastic framework to be used in remainder of the paper.

Stochastic Framework
In this paper we will focus on the stochastic framework introduced in [25] referred to as Stochastic Discrete Dynamical Systems (SDDS).This framework is a natural extension of Boolean networks and is an appropriate setup to model the effect of intrinsic noise on network dynamics.Consider the discrete variables x 1 , . . ., x n that can take values in finite sets S 1 , . . ., S n , respectively.Let S = S 1 × • • • × S n be the Cartesian product.A SDDS in the variables x 1 , . . ., x n is a collection of n triplets where • f i : S → S i is the update function for x i , for all i = 1, . . ., n.
• p ↑ i is the activation propensity.
• p ↓ i is the degradation propensity.
. These are the parameters of interest in this paper.
The stochasticity originates from the propensity parameters p ↑ k and p ↓ k , which should be interpreted as follows: If there would be an activation of x k at the next time step, i.e., if s 1 , s 2 ∈ S k with s 1 < s 2 and x k (t) = s 1 , and f k (x 1 (t), . . ., x n (t)) = s 2 , then x k (t + 1) = s 2 with probability p ↑ k .The degradation probability p ↓ k is defined similarly.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 Let x, y ∈ S. The transition from x to y is given by Notice that P rob(x i → y i ) = 0 for all y i / ∈ {x i , f i (x)}.Then the transition matrix is given by The dynamics of SDDS depends on the transition probabilities a xy , which depend on the propensity values and the update functions.Online software to test examples is available at http://adam.plantsimlab.org/(choose Discrete Dynamical Systems (SDDS) in the model type).
In Markov chain notation, the transition probability a xy = p(X t = x|X t−1 = y) represents the probability of being in state x at time t given that system was in state y at time t − 1.If π t = p(X t = x) represents the probability of being in state x at time t, then we will assume that π is a row vector containing the probabilities of being in state x at time t for all x ∈ S. If π 0 is the initial distribution at time t = 0, then at time t = 1, If we iterate Equation 3 and if we get to the point where then we will say that the Markov chain has reached a stationary distribution and that π is the stationary distribution.

Methods
In this section we describe a method for estimating the propensity parameters for SDDS.The approach is based on adding noise to the system using the Google PageRank [4,16,23] strategy to make the system ergodic and thus guaranteeing the existence of a stationary distribution and then with the use of a genetic algorithm the propensity parameters are estimated.To guarantee the existence of a stationary distribution we use a special case of the Perron-Frobenius Theorem.
• The vector π is the unique probability vector which is an eigenvector of A associated with the eigenvalue 1.
A proof of Theorem 2.1 can be found in Chapter 10 of [16].Theorem 2.1 ensures a unique stationary distribution π provided that we have a regular transition matrix, that is, if some power A k contains only strictly positive entries.However, the transition matrix A of SDDS given in Equation 2might not be regular.In the following subsection, we use a similar approach to the Google's PageRank algorithm to add noise to the system to obtain a new transition matrix that is regular.

PageRank Algorithm
For simplicity, consider a SDDS, F = {f i , p ↑ i , p ↓ i } n i=1 where f i : S → S i , S = K n , and |K| = p.Then its transition matrix A given in Equation 2 is a p n × p n matrix.To introduce noise into the system we consider the Google Matrix where g is a constant number in the interval [0, 1] and K is a p n × p n matrix all of whose columns are the vector (1/p n , . . ., 1/p n ).The matrix G in Equation 5 is a regular matrix and then we can use Theorem 2.1 to get a stationary distribution for G, This stationary distribution reflects the dynamics of the SDDS The importance of a state x ∈ S can be measured by the size of the corresponding entry π x in the stationary distribution of Equation 6.For instance, for ranking the importance of the states in a Markov chain one can use the size of the corresponding entries in the stationary distribution.We will refer to this entry π x as the PageRank score of x.

Genetic Algorithm
The entries of the stationary distribution π in Equation 6 can also be interpreted as occupation times for each state.Thus it gives the probability of being at a certain state.Now suppose that we start with a desired stationary distribution π * = (π * 1 , . . ., π * p n ).We have developed a genetic algorithm that initializes a population of random propensity matrices and searches for a propensity matrix prop * such that its stationary distribution π = (π 1 , . . ., π p n ) gets closer to the desired stationary distribution π * .That is, we search for propensity matrices such that the distance between π and π * is minimized, min The pseudocode of this genetic algorithm is given in Algorithm 1 and it has been implemented in Octave/Matlab and our code can be downloaded from http://www.ms.uky.edu/∼dmu228/GeneticAlg/Code.html.
Algorithm 1 Genetic Algorithm with PageRank.

Estimating the stationary distribution
The genetic algorithm, Algorithm 1, uses the exact stationary distribution through PageRank (see Equation 6) which is computationally expensive for larger models.Here we present an efficient algorithm for estimating the stationary distribution based on a random walk.The expensive part of Algorithm 1 is the calculation of the stationary distribution π in Equation 6.We have implemented an algorithm for estimating the stationary distribution by doing a random walk using SDDS as a generative model; see Algorithm 2. The idea behind Algorithm 2 is to use SDDS for simulating from an initial state according to the transition probabilities given in Equation 5.That is, we initialize the simulation at an initial state x ∈ S and then with probability g we move to another state y ∈ S according to a xy (see Equation 1) and with probability 1 − g we jump to a random node.We repeat this process for a given number of iterations.At the end of a maximum number of iterations, we count how often we visited each state and the normalized frequencies will be the approximated stationary distribution.
To make the genetic algorithm more efficient, we used the estimated stationary distribution described in the previous paragraph.Thus, within the fitness function of Algorithm 1, we use the estimated stationary distribution to assess the fitness of the generated propensity matrices.The pseudocode for this new algorithm is the same as Algorithm 1, the only change is in the fitness function.This version of the algorithm has also been implemented in Octave/Matlab and our code can be found in http://www.ms.uky.edu/∼dmu228/GeneticAlg/Code.html.

Results
We test our methods using two published models that are appropriate for changing the stationary distribution under the choice of different propensity parameters.The first model is a Boolean network while the second is a multistate model.In both models bistability has been observed but the basin size of one of the attractors under the synchronous update is much larger than the basin of the other attractor, and thus the stationary distribution will be more concentrated in one of the attractors.We will use our methods to change the stationary distribution in favor of the attractor with a smaller basin.
Example 2.2.Lac-operon network.The lac-operon in E. coli [12] is one of the best studied gene regulatory networks.This system is responsible for Algorithm 2 Estimate Stationary Distribution.

6:
for i=1,. . ., N umIter do if rand < g then 8: y = random state between 1 and p n .else  return π the metabolism of lactose in the absence of glucose.This system exhibits bistability in the sense that the operon can be either ON or OFF, depending on the presence of the preferred energy source: glucose.A Boolean network for this system has been developed in [37].This model considers the following 10 components x 1 = M: lac mRNA, x 2 = P: lac permease, x 3 = B: lacβ-galactosidase, x 4 = C: CAP, x 5 = R: repressor, x 6 = Rm: repressor at medium concentration, x 7 = A: allolactose, x 8 = Am: allolactose at medium concentration, x 9 = L: lactose, x 10 = Lm: lactose at medium concentration, and the Boolean rules are given by where G e , L em , and L e are parameters.G e and L e indicate the concentration of extracellular glucose and lactose, respectively.The parameter L em indicates the medium concentration of extracellular lactose.For medium extracellular lactose, that is when L em = 1 and L e = 0, this system has two fixed points: s 1 = (0, 0, 0, 1, 1, 1, 0, 0, 0, 0) and s 2 = (1, 1, 1, 1, 0, 0, 0, 1, 0, 1) that represent the state of the operon being OFF and ON, respectively.To test our method we calculated the stationary distribution of the system using Equation 6with g = 0.9 in Equation 5. First we used the propensity values given in Equation 10where all propensities are fixed to 0.9.This choice of parameters approximates the synchronous dynamics in the sense that each function has a 90% change of being used during the simulations and 10% chance of maintaining its current value.Under this selection of parameters, the fixed point s 1 has a PageRank score of 0.3346 while the other fixed point s 2 has a score of 0.0463, see Figure 1a and Table 1.Then we have applied our genetic algorithm to search for parameters that can increase the PageRank score of s 2 and decrease the score of s 1 .After doing so, we found the propensity parameters given in Equation 11.With this new set of parameters, the fixed point s 1 has a score of 0.0199 while s 2 has a score of 0.5485, see Figure 1b and Table 1.To appreciate the impact of the change in propensity parameters, we have plotted the state space of the system with both propensity matrices in Figure 2. The edges in blue in Figure 2 represent the most likely trajectory.Notice that in Figure 2b the trajectories are leading towards s 2 and the size of the labels of the nodes were scaled according to their PageRank score, see Table 1.
(b) Scores with propensities in Equation 11.
Figure 1: PageRank scores before and after the genetic algorithm.In each panel, the x-axis shows the PageRank scores while the y-axis shows the frequencies of states with the given scores in the x-axis (the exact scores for the states of interest are given in Table 1).Left panel shows the state space where all the propensities are equal to 0.9 while the right panel shows the state space where the propensity parameters where estimated using the genetic algorithm.

Propensities
Attractor Score In Equation 100001110000 0.3346 (all fixed to 0.9) 1111000101 0.0463 In Equation 110001110000 0.0199 (genetic algorithm) 1111000101 0.5485 Table 1: PageRank scores for the states of the attractors of the system.The order of variables in each vector state is M, P, B, C, R, R m , A, A m , L, L m .
Example 2.3.Phage lambda infection.The outcome of phage lambda infection is another system that has been widely studied over the last decades [13,26,32,33,38].One of the earliest models that has been developed for this  Left panel shows the state space where all the propensities are equal to 0.9 while the right panel shows the state space with the estimated propensity parameters using the genetic algorithm.The edges in blue represent the most likely trajectory.The size of the labels of the nodes were scaled according to their PageRank score.
system is the logical model by Thieffry and Thomas [33].The regulatory genes considered in is this model are: CI, CRO, CII, and N. Experimental reports [15,28,32,33] have shown that, if the gene CI is fully expressed, all other genes are OFF.In the absence of CRO protein, CI is fully expressed (even in the absence of N and CII).CI is fully repressed provided that CRO is active and CII is absent.
This network is a bistable switch between lysis and lysogeny.Lysis is the state where the phage will be replicated, killing the host.Otherwise, the network will transition to a state called lysogeny where the phage will incorporate its DNA into the bacterium and become dormant.These cell fate differences have been attributed to the spontaneous changes in the timing of individual biochemical reaction events [20,33].
In the model of Thieffry and Thomas [33], the first variable, CI, has three levels {0, 1, 2}, the second variable, CRO, has four levels {0, 1, 2, 3}, and the third and fourth variables, CII and N , are Boolean.Since the nodes of this model have different number of states, in order to apply our methods, we have extended the model so that all nodes have the same number of states.We have used the method given in [36] to extend the number of states such that all nodes have 5 states (the method for extending requires a prime number for number of states so we have chosen 5 states).The method for extending the number of states preserves the original attractors.The update rules for this model are available with our code that is freely available.The extended model has a steady state, 2000, and a 2-cycle involving 0200 and 0300.The steady state 2000 represents lysogeny where CI is fully expressed while the other genes are OFF.The cycle between 0200 and 0300 represents lysis where CRO is active and other genes are repressed.
To test our method we calculated the stationary distribution of the system using Equation 6with g = 0.9 in Equation 5. First we used the propensity values given in Equation 12where all propensities are fixed to 0.9.This choice of parameters approximates the synchronous dynamics in the sense that each function has a 90% change of being used during the simulations and 10% chance of maintaining its current value.Under this selection of parameters, the fixed point 2000 has a PageRank score of 0.2772 while the states of the cycle 0200 and 0300 have scores of 0.2185 and 0.2108, respectively.Notice that this cycle attractor will have an overall score of 0.4293, see Figure 3a and Table 2. Then we have applied our genetic algorithm to search for parameters that can increase the PageRank score of the fixed point 2000 and found the propensity parameters given in Equation 13.With this new set of parameters, the fixed point 2000 has a score of 0.6040 while the states of the cycle 0200 and 0300 have scores of 0.0716 and 0.00016, respectively, see Figure 3b and Table 2. To appreciate the impact of the change in propensity parameters, we have plotted the state space of the system with both propensity matrices in Figure 4.The edges in blue in Figure 4 represent the most likely trajectory.Notice that in Figure 4b the trajectories are leading towards 2000 and the size of the labels of the nodes were scaled according to their PageRank score, see Table 2.
CI CRO CII N p ↑ i .9 .9 .9.9p ↓ i .9 .9 .9.92).Left panel shows the PageRank scores where all the propensities were equal to 0.9 while the right panel shows the scores where the propensity parameters where estimated using the genetic algorithm.The score for the state 2000 is 0.6040.

Propensities
Attractor Score In

Discussion
Parameter estimation for stochastic models of biological networks is in general a very hard problem.This paper focuses on a class of stochastic discrete models, which is an extension of Boolean networks.The methods presented here use a well stablished approach for introducing noise into a system in a way that ergodicity and thus the existence of a unique stationary distribution is guaranteed.Then a genetic algorithm for calculating a set of parameters able to approximate a desired stationary distribution was developed.Also, techniques for approximating the stationary distribution at each iteration of the genetic algorithm that make the search process more efficient was applied.
One shortcoming of the method is that it is a stochastic method.That is, each time we the algorithm we get a different result.An exhaustive investigation about the variance of the results is still missing.For the examples that we presented in the results section, the output of the algorithm might vary in about 20% of the reported propensities.
In parameter estimation, sometimes, it is useful to identify a smaller set of key parameters to estimate.This problem is out of the scope of the paper.However, for Boolean network models, one way to address this problem could be by using the different network reduction algorithms, for instance see [29,35], to identify a smaller "core" network that preserves the important features of the dynamics of the original network.And then one could apply the parameter estimation techniques that are described in this paper.This type of approach could be especially useful if dealing with very large networks where running the genetic algorithm is computationally expensive.

Conclusions
In this paper we present an efficient method for estimating the parameters of a stochastic framework.The modeling framework is an extension of Boolean networks that uses propensity parameters for activation and inhibition.Parameter estimation techniques are needed whenever one needs to tune the propensity parameters of the stochastic system to reproduce a desired stationary distribution.For instance, if dealing with a bistable system and if it is desired to have the stationary distribution that have the PageRank score concentrated in one of the attractors of the system, then one needs to estimate the propensity parameters that represent such a desired distribution.Parameter estimation methods for this purpose were not available.In this paper we present a method for estimating propensity parameters given a desired stationary distribution for the system.We tested the method in one Boolean network with 10 nodes (where the size of the state space is 2 10 = 1024) and a multistate network with 4 nodes where each node has 5 states (where the size of the state space is 4 5 = 625).For each system, we were able to redirect the system towards the attractor with the smaller PageRank score.The method is efficient and for the examples we have shown it can be run in few seconds in a laptop computer.Our code is available at http://www.ms.uky.edu/∼dmu228/GeneticAlg/Code.html.

Figure 2 :
Figure 2: State space comparison before and after the genetic algorithm.Left panel shows the state space where all the propensities are equal to 0.9 while the right panel shows the state space with the estimated propensity parameters using the genetic algorithm.The edges in blue represent the most likely trajectory.The size of the labels of the nodes were scaled according to their PageRank score.
Scores with propensities in Equation12.(b)Scores with propensities in Equation13.

Figure 3 :
Figure 3: PageRank scores before and after the genetic algorithm.In each panel, the x-axis shows the PageRank scores while the y-axis shows the frequencies of states with the given scores in the x-axis (the exact scores for the states of interest are given in Table2).Left panel shows the PageRank scores where all the propensities were equal to 0.9 while the right panel shows the scores where the propensity parameters where estimated using the genetic algorithm.The score for the state 2000 is 0.6040.

Figure 4 :
Figure 4: State space comparison before and after the genetic algorithm.Left panel shows the state space where all the propensities are equal to 0.9 while the right panel shows the state space with the estimated propensity parameters using the genetic algorithm.The edges in blue represent the most likely trajectory.The size of the labels of the node were scaled according to their PageRank score.

Table 2 :
PageRank scores for the states of the attractors of the system.The order of variables in each vector state is CI, CRO, CII, N .