Gene perturbation and intervention in context-sensitive stochastic Boolean networks
- Peican Zhu^{1},
- Jinghang Liang^{1} and
- Jie Han^{1}Email author
DOI: 10.1186/1752-0509-8-60
© Zhu et al.; licensee BioMed Central Ltd. 2014
Received: 11 February 2013
Accepted: 22 April 2014
Published: 21 May 2014
Abstract
Background
In a gene regulatory network (GRN), gene expressions are affected by noise, and stochastic fluctuations exist in the interactions among genes. These stochastic interactions are context dependent, thus it becomes important to consider noise in a context-sensitive manner in a network model. As a logical model, context-sensitive probabilistic Boolean networks (CSPBNs) account for molecular and genetic noise in the temporal context of gene functions. In a CSPBN with n genes and k contexts, however, a computational complexity of O(nk ^{2}2^{2n }) (or O(nk 2^{ n })) is required for an accurate (or approximate) computation of the state transition matrix (STM) of the size (2^{ n } ∙ k) × (2^{ n } ∙ k) (or 2^{ n } × 2^{ n }). The evaluation of a steady state distribution (SSD) is more challenging. Recently, stochastic Boolean networks (SBNs) have been proposed as an efficient implementation of an instantaneous PBN.
Results
The notion of stochastic Boolean networks (SBNs) is extended for the general model of PBNs, i.e., CSPBNs. This yields a novel structure of context-sensitive SBNs (CSSBNs) for modeling the stochasticity in a GRN. A CSSBN enables an efficient simulation of a CSPBN with a complexity of O(nLk 2^{ n }) for computing the state transition matrix, where L is a factor related to the required sequence length in CSSBN for achieving a desired accuracy. A time-frame expanded CSSBN can further efficiently simulate the stationary behavior of a CSPBN and allow for a tunable tradeoff between accuracy and efficiency. The CSSBN approach is more efficient than an analytical method and more accurate than an approximate analysis.
Conclusions
Context-sensitive stochastic Boolean networks (CSSBNs) are proposed as an efficient approach to modeling the effects of gene perturbation and intervention in gene regulatory networks. A CSSBN analysis provides biologically meaningful insights into the oscillatory dynamics of the p53-Mdm2 network in a context-switching environment. It is shown that random gene perturbation has a greater effect on the final distribution of the steady state of a network compared to context switching activities. The CSSBN approach can further predict the steady state distribution of a glioma network under gene intervention. Ultimately, this will help drug discovery and develop effective drug intervention strategies.
Keywords
Gene regulatory networks Boolean networks Stochastic Boolean networks Context dependent Gene perturbation Intervention Context switch Steady state distribution p53 network glioma networkBackground
Diverse biological functions are regulated through the interactions among genes, proteins and other molecules in a cell. Gene expressions are however affected by the intrinsic and extrinsic noise in a gene network [1]. A major source of the noise is the stochastic fluctuations in gene regulatory interactions [2]. The genetic interactions are also context dependent, that is, certain regulatory functions are active in some cellular states, but inactive in others [3]. This indicates the necessity to consider noise in a context-sensitive manner in the study of gene regulatory networks (GRNs).
Various methods have been proposed to model GRNs; these include logical models [4], continuous models using differential equations [5, 6] and stochastic models at the single-molecule level [7, 8]. As a classic logical model, Boolean networks (BNs) have been widely used to qualitatively model the interactions among genes [4, 9–11]. Probabilistic Boolean networks (PBNs) have been proposed to consider noise in a BN model [12–14]. In a PBN, the next state of a gene is determined by its current state and a Boolean function. If the Boolean function is randomly selected, a PBN is referred to as an instantaneous PBN [12]. As a general model, a context-sensitive PBN (CSPBN) considers the feature of context dependence in a biological network [15]. In a CSPBN, a context is a combination of Boolean functions and each function determines the next state of a gene. A context remains unchanged until a switching occurs. This switching of contexts, possibly caused by external stimuli, is considered to occur randomly in a network.
The study of PBNs has focused on the analysis of steady state distributions (SSDs) under gene perturbation and intervention. A Markov Chain Monte Carlo method is used in [16] to analyze the long run behavior of a PBN; however, this method generally requires a large number of simulations to reach a steady state, due to the slow convergence typically encountered in a Monte Carlo method [17]. An analysis is performed in [18] for finding the SSD of a PBN through the computation of the state transition matrix (STM). However, the application of an analytical approach is generally limited to small networks due to the exponential increase in the size of an STM with gene numbers. The analysis of CSPBNs presents even a greater challenge due to its significantly increased computational complexity. Analytical expressions have been derived for analyzing CSPBNs [19]; however, this method is only applicable to the steady state analysis of a network with small perturbation and switching probabilities. The method in [20] ignores some BNs with very small probabilities for reducing the size of the STM and thus provides a more efficient but approximate solution for computing the SSD of a CSPBN. In [15, 21, 22], gene intervention is investigated for avoiding undesirable states associated with certain diseases (such as cancer). Due to external stimuli, the STM is changed by external control variables, so desirable states can be obtained with larger probabilities in the SSD. In a context-sensitive network with n genes and k contexts, however, a (2^{ n } ∙ k) × (2^{ n } ∙ k) [23] (or 2^{ n } × 2^{ n }[20]) matrix is required for an accurate (or approximate) analysis of the SSD; this results in a computational complexity of O(nk ^{2}2^{2n}) (or O(nk 2^{ n })) for an accurate (or approximate) computation of the STM. Hence, the application of current gene network analysis is limited to those of less than a dozen of genes.
Recently, stochastic Boolean networks (SBNs) have been proposed for an efficient computation of the STM and SSD of an instantaneous PBN [24]. The SBN approach can recover biologically-proven regulatory behaviors, such as the oscillatory dynamics of a simplified p53-Mdm2 network [25] and the dynamic attractors in a T cell immune response network [26]. To further exploit the simplicity of logical models, stochastic multiple-valued networks (SMNs) have been utilized to investigate the dynamics of gene networks with multiple-valued gene states [27]. Asynchronous SBNs have also been developed to investigate the asynchronous state-update behavior of genes [28]. In this paper, the notion of SBNs is extended for the general model of PBNs, i.e., the CSPBNs. This presents a novel structure of context-sensitive SBNs (CSSBNs) for modeling the stochasticity in a context-dependent GRN. In particular, gene perturbation and intervention are considered in CSSBNs. Through an efficient simulation of a CSSBN, the computational complexity in the evaluation of a CSPBN is reduced from O(nk ^{2}2^{2n }) to O(nLk 2^{ n }) for computing the STM, where L is a factor related to the required sequence length in CSSBN for achieving a desired accuracy. The use of non-Bernoulli sequences of random permutations of fixed number of 1’s and 0’s further increases computational efficiency and allows for a tunable tradeoff between accuracy and efficiency. A time-frame expanded CSSBN can further simulate the stationary behavior of a CSPBN and produce results that are more accurate than an approximate analysis, thus making the CSSBN useful in modeling complex context-sensitive GRNs. The CSSBN models are applied to the study of a simplified p53-Mdm2 network and the results using an SSD control policy are reported on the effect of external gene intervention in a glioma network [18, 20].
Methods
Context-Sensitive Probabilistic Boolean Networks (CSPBNs)
For a network of n genes, a probabilistic Boolean network (PBN) is defined by G(V, F), where V = {x _{1}, x _{2}, …, x _{ n }}, a set of binary-valued nodes, F = (F _{1}, F _{2}, …, F _{ n }), a list of sets of Boolean functions: ${\mathit{F}}_{\mathit{i}}=\left\{{\mathit{f}}_{1}^{\left(\mathit{i}\right)},{\mathit{f}}_{2}^{\left(\mathit{i}\right)},\dots ,{\mathit{f}}_{\mathit{l}\left(\mathit{i}\right)}^{\left(\mathit{i}\right)}\right\}$ and l(i) is the number of possible functions for gene i, i = 1, 2, …, n[12–14]. A node x _{ i } represents the state of gene i; x _{ i } = 1 (or 0) indicates that gene i is (or not) expressed. The set F _{ i } contains the rules that determine the next state of gene i. Each ${\mathit{f}}_{\mathit{j}\left(\mathit{i}\right)}^{\left(\mathit{i}\right)}:{\left\{0,1\right\}}^{\mathit{n}}\to \left\{0,1\right\}$ for 1 ≤ j(i) ≤ l(i) is a mapping or predictor function determining the state of gene i. In a context-sensitive PBN (CSPBN) with k contexts, a network function is defined as ${\mathit{f}}_{\mathit{j}}=\left({\mathit{f}}_{\mathit{j}}^{\left(1\right)},{\mathit{f}}_{\mathit{j}}^{\left(2\right)},\dots ,{\mathit{f}}_{\mathit{j}}^{\left(\mathit{n}\right)}\right)$, where ${\mathit{f}}_{\mathit{j}}^{\left(\mathit{i}\right)}:{\left\{0,1\right\}}^{\mathit{n}}\to \left\{0,1\right\}$ is the predictor function for gene i, i = 1, 2, …, n, in context j (j = 1, 2, …, k) [15, 29].
Due to the stochasticity in genetic networks, the next state of gene i is determined by all the l(i) functions in F _{ i }, i.e., by ${\mathit{f}}_{1}^{\left(\mathit{i}\right)},{\mathit{f}}_{2}^{\left(\mathit{i}\right)},\dots ,{\mathit{f}}_{\mathit{l}\left(\mathit{i}\right)}^{\left(\mathit{i}\right)}$, with probabilities ${\mathit{c}}_{1}^{\left(\mathit{i}\right)},{\mathit{c}}_{2}^{\left(\mathit{i}\right)},\dots ,{\mathit{c}}_{\mathit{l}\left(\mathit{i}\right)}^{\left(\mathit{i}\right)}$. For an instantaneous PBN of n genes, there are a total number of $\prod _{\mathit{i}=1}^{\mathit{n}}}\mathit{l}\left(\mathit{i}\right)$ possible Boolean networks (BNs), each of which is a possible realization of the genetic network. Each BN can be considered as the network function for a context [20]. In the j th context, assume the network function is given by ${\mathit{f}}_{\mathit{j}}=\left({\mathit{f}}_{\mathit{j}\left(1\right)}^{\left(1\right)},{\mathit{f}}_{\mathit{j}\left(2\right)}^{\left(2\right)},\cdots ,{\mathit{f}}_{\mathit{j}\left(\mathit{n}\right)}^{\left(\mathit{n}\right)}\right)\phantom{\rule{0.12em}{0ex}}$, where each ${\mathit{f}}_{\mathit{j}\left(\mathit{i}\right)}^{\left(\mathit{i}\right)}:{\left\{0,1\right\}}^{\mathit{n}}\to \left\{0,1\right\}$ for 1 ≤ j(i) ≤ l(i) is a mapping or predictor function determining the state of gene i. The probability that the j th context is selected, is obtained as ${\mathit{C}}_{\mathit{j}}={\displaystyle \prod _{\mathit{i}=1}^{\mathit{n}}}{\mathit{c}}_{\mathit{j}\left(\mathit{i}\right)}^{\left(\mathit{i}\right)}$ for j = 1, 2, …, k, where k is the number of contexts in a context-sensitive network [20]. The state of a gene is updated in the selected context; thus the next state depends on both the present state and the selected context.
In a context-sensitive network, a context may remain for certain time until a random event occurs. A context switching usually occurs with probability q. For q = 1, the CSPBN becomes an instantaneous PBN. For q < 1, if a new context is to be selected, it is randomly chosen from the set of network functions: {f _{1}, f _{2}, …, f _{ k }}, with a set of context selection probabilities: {C _{1}, C _{2}, …, C _{ k }} [29]. If noise is considered in a CSPBN, it is often referred to as a perturbation, by which a gene flips its state with a probability p (p ≠ 0). It has been shown that a PBN with perturbation is an ergodic Markov chain in that all the states are connected in the PBN [13]. The transition probability for any two states is determined by the values of p and q pairs. Following [19], one of four mutually exclusive events occurs at time t in a CSPBN with perturbation:
∅_{1}: The predictor functions in the currently selected context are applied to update the gene expressions and this context remains for the next transition.
∅_{2}: The predictor functions in the currently selected context are applied to update the states of the genes and then a new context is selected for the next transition.
∅_{3}: A random perturbation occurs and the currently selected context remains for the next transition.
∅_{4}: A random perturbation occurs and a new context is selected for the next transition.
The effect of a switching order, i.e., whether a network switches its context before or after its state transition, is considered in [23], whereas in this paper, we focus on the transition rules, i.e., the network function is applied first and then the context switches.
A gene activity profile (GAP) is defined as a vector for describing the state of a network at time t, x(t) = (x _{1}(t), x _{2}(t), …, x _{ n }(t)), where x _{ i }(t) ∈ {0, 1} for i = 1, 2, …, n. The state of a CSPBN can be represented as a combination of a context and a GAP, i.e., S = (context j, GAP i). By this definition, the number of states in a CSPBN increases from 2^{ n } to 2^{ n } k for a network with n genes and k contexts. A GAP can also be represented by its decimal index, i.e., $\mathit{d}={\displaystyle \sum _{\mathit{j}=1}^{\mathit{n}}}{2}^{\mathit{n}-\mathit{j}}{\mathit{x}}_{\mathit{j}}\left(\mathit{t}\right)+1$. A state of a CSPBN is then given by S = (f _{ i }, d) for i ∈ {1, 2, …, k} and d ∈ {1, 2, …, 2^{ n }}, where d is the decimal index of a GAP. For convenience, a state (f _{ i }, d) is referred to as (i, d) in the following analysis.
and h is the Hamming distance between two GAPs with decimal indices x _{1} and x _{2}, x _{1}, x _{2} ∈ {1, 2, …, 2^{ n }}. This distance indicates the number of genes with different expressions in the two GAPs x _{1} and x _{2}.
where P _{ i } and ${\tilde{\mathit{P}}}_{\mathit{i}}$ denote the STMs for the Boolean network i, i ∈ {1, 2, …, k}, without and with perturbation respectively.
Context-Sensitive Stochastic Boolean Networks (CSSBNs)
A CSSBN without perturbation
s and r denote the s th and r th contexts, y and x represent two gene activity profiles (GAPs, in decimal indices), and C _{ s } and C _{ r } indicate the probability of selecting the s th and r th contexts respectively.
If a switch does occur, the context selection process is implemented by another multiplexer for choosing one from the k contexts, according to predefined selection probabilities. As a network function is a combination of each gene’s predictor function, a combination of the m control sequences of S _{1} ~ S _{ m } is used to encode the predefined selection probabilities; this in turn determines the selection probability of each new context. For a CSPBN with n genes, the CSSBN needs to be run for each of the 2^{ n } input states and sequences need to be generated for the m control signals of the multiplexer.
For a switching probability of q in the proposed CSSBN, the currently selected context i remains at time t with a probability of q and switches to one of the k contexts with a probability of 1 - q. If the network transitions from a GAP y to x in the s th context (i.e., f _{ s,y,x } = 1), then context s will remain with probability 1 - q + qC _{ s } at the next time step. Otherwise, the network moves into a new context r with probability qC _{ r }. From this analysis, it can be seen that the CSSBN in Figure 2 computes the transition probability from state (s, y) to (r, x) as given by (7). This indicates that the proposed CSSBN model accurately implements the function of a CSPBN.
A CSSBN with perturbation
which is equivalent to (9). This indicates that a CSPBN with perturbation can be accurately implemented by a CSSBN with perturbation.
Intervention in a CSSBN
In contrast to perturbation, gene intervention refers to the process of deliberately changing the states of some genes to guide a network into a desired state [12]. External stimuli are applied to a network to avoid undesirable states that might be associated with certain diseases. For an effective intervention, control policies are developed for different intervention strategies; as several genes may affect the state of the target gene, a single gene that is most influential on the network state is usually identified as the control gene, due to its simplistic biological implications [15]. The STM is then changed by an external intervention [15, 21, 22].
In a CSPBN, the context information is usually hidden and thus difficult to obtain in practice. However, gene expressions can readily be observed, so a GAP within all possible contexts is considered to be undesirable if the GAP is an undesirable state of the network [36]. Hence, gene intervention refers to changing the state of a network as represented by a GAP.
For an n-gene network, a vector s = (s _{1}, s _{2}, s _{3}, …, s _{ n }) with s _{ i } ∈ {0, 1} for any i ∈ {1, 2, …, n}, is defined as the control gene vector, that is, if s _{ i } = 1, then gene i is selected as a control gene [37]. An intervention vector is defined as $\mathit{u}=\left({\mathit{u}}_{1},{\mathit{u}}_{2},{\mathit{u}}_{3},\dots ,{\mathit{u}}_{{2}^{\mathit{n}}}\right)$, u _{ j } ∈ {0, 1} for any j ∈ {1, 2, …, 2^{ n }}, where u _{ j } = 1 (or 0) indicates a flipping (or remaining) of the state of a control gene at GAP j. If gene i is selected to be a control gene, for example, then s _{ i } = 1. The expression level of the control gene i is then determined by its current state and the status of u. If u _{ j } = 1, i.e., s _{ i } u _{ j } = 1, the state of control gene i is flipped by the external intervention when GAP j emerges; otherwise, the state of control gene i remains unchanged and thus the current state is preserved in GAP j[38, 39]. An intervention vector can be obtained by various methods; in this paper, a control policy using the steady state distribution (SSD) [40] is used to obtain u.
where s _{ g } is a vector of n bits with 0 at every bit except the g th bit. 1(∙) is an indicator function: $1\left(\mathit{u}\left({\widehat{\mathit{x}}}_{\mathit{t}}\right)=1\right)=1$ if $\mathit{u}\left({\widehat{\mathit{x}}}_{\mathit{t}}\right)=1$ and $1\left(\mathit{u}\left({\widehat{\mathit{x}}}_{\mathit{t}}\right)=1\right)=0$ otherwise. When $\mathit{u}\left({\widehat{\mathit{x}}}_{\mathit{t}}\right)=0$, the GAP ${\widehat{\mathit{x}}}_{\mathit{t}}$ is maintained; otherwise, the GAP ${\widehat{\mathit{x}}}_{\mathit{t}}$ is an undesirable state, for which the control bit is to be toggled, thereby improving the occurring probability of a desirable state.
which is equivalent to (11). This indicates that the CSSBN model in Figure 4 accurately implements the process of external gene intervention.
State transition matrix (STM) and steady state analysis
Minimum sequence length and average run time in computing the state transition matrix (STM) for context-sensitive stochastic Boolean networks (CSSBNs)
n | k | Num | CSSBN (Norm 2 = 0.05) | CSSBN (Norm 2 = 0.03) | CSPBN[23] | |||||
---|---|---|---|---|---|---|---|---|---|---|
L(bits) | Average time(s) | Standard deviation | L(bits) | Average time(s) | Standard deviation | Average time(s) | Standard deviation | |||
2 | 1 | 4 | 40 | 0.00148 | 0.00048 | 100 | 0.00181 | 0.00048 | 0.00049 | 0.00046 |
2 | 4 | 16 | 900 | 0.00553 | 0.00167 | 2000 | 0.01116 | 0.00036 | 0.01248 | 0.00251 |
2 | 16 | 64 | 1000 | 0.64383 | 0.41941 | 6000 | 1.02666 | 0.10779 | 0.09593 | 0.00901 |
3 | 1 | 8 | 140 | 0.00446 | 0.00132 | 340 | 0.00812 | 0.00094 | 0.00157 | 0.00045 |
4 | 1 | 16 | 250 | 0.01438 | 0.00287 | 750 | 0.03685 | 0.00268 | 0.00443 | 0.00027 |
4 | 16 | 256 | 8000 | 1.43115 | 0.01844 | 20000 | 3.29188 | 0.03830 | 1.66014 | 0.10349 |
5 | 1 | 32 | 420 | 0.04602 | 0.00261 | 1200 | 0.12341 | 0.00314 | 0.01862 | 0.00214 |
5 | 32 | 1024 | 8000 | 10.7359 | 0.01350 | 25000 | 35.5412 | 0.04335 | 24.0659 | 0.08390 |
6 | 1 | 64 | 600 | 0.09993 | 0.00458 | 1700 | 0.33797 | 0.00197 | 0.06662 | 0.00705 |
6 | 64 | 4096 | 24000 | 74.5058 | 1.20649 | 60000 | 195.793 | 4.86945 | 531.068 | 6.36186 |
7 | 1 | 32 | 1000 | 0.52519 | 0.02762 | 2800 | 1.34503 | 0.06377 | 0.25873 | 0.00719 |
7 | 16 | 2048 | 8000 | 18.8938 | 0.42939 | 25000 | 50.0857 | 0.27419 | 221.922 | 1.65582 |
8 | 1 | 256 | 1400 | 1.45711 | 0.02798 | 3500 | 3.67661 | 0.16084 | 1.05147 | 0.01302 |
8 | 4 | 1024 | 4000 | 5.50577 | 0.05167 | 12500 | 15.6728 | 0.11539 | 37.4102 | 2.20907 |
9 | 1 | 512 | 2600 | 5.36480 | 0.07636 | 8000 | 16.2693 | 0.14106 | 4.06749 | 0.06899 |
9 | 4 | 2048 | 8000 | 23.1989 | 0.40195 | 18000 | 45.9247 | 0.78522 | 72.4459 | 2.62472 |
10 | 1 | 1024 | 4000 | 18.4761 | 0.12540 | 11000 | 49.2457 | 0.36526 | 17.5134 | 0.12540 |
10 | 4 | 4096 | 10000 | 59.9416 | 0.41164 | 30000 | 162.686 | 1.13288 | 619.311 | 13.4098 |
In general, a (2^{ n } ∙ k) × (2^{ n } ∙ k) STM is required for an accurate CSPBN analysis [23] while a 2^{ n } × 2^{ n } STM is needed by an approximate method [20]. As the number of genes increases in a network, a matrix-based analysis becomes infeasible due to a significant increase in the required computational resources. The analysis of the SSD is even more challenging for a CSPBN. For a CSSBN, however, the STM can be accurately and efficiently computed. Furthermore, the SSD can be evaluated through an iterative simulation in the temporal domain (or the so-called time-frame expansion technique [24]). By this technique, an iterative structure of the CSSBN is used to simulate the temporal evolution of a GRN. The required number of iterations is determined by the number of state transitions before reaching a steady state.
where ϵ indicates a threshold for determining whether the steady state has been reached or not. If (14) is met, a network is considered to have reached a steady state, i.e., x ^{(∞)} = x ^{(m)}; thus x ^{(m)} is considered as the stationary distribution. Alternatively, a reasonable number of network transitions (e.g. a few hundreds) can be run to obtain the SSD in practice. As shown in the Results and Discussion section, a time-frame extended simulation of CSSBNs is more efficient than an analytical approach while producing more accurate results than an approximate analysis. See Additional file 2: The pseudocode for computing the SSD of a CSSBN.
In a time-frame expanded CSSBN, random binary bit streams are generated for predictor function selection, gene perturbation and context switching probabilities. As in [24], non-Bernoulli sequences of random permutations of fixed numbers of 1’s and 0’s are used for encoding initial probabilities. These sequences then propagate through the CSSBNs and the statistics encoded in the output sequences are used to obtain the SSD. Compared to the use of Bernoulli sequences, as shown in [35], the use of the non-Bernoulli sequences produces more accurate results for a given sequence length or requires a shorter sequence length for a desired evaluation accuracy. As an efficient alternative to an STM-based analysis, the stochastic simulation of CSSBNs provides flexibility in achieving a tunable accuracy-efficiency tradeoff by using stochastic sequences of different lengths. Hence, the proposed CSSBN approach is potentially useful in the analysis of large GRNs.
Results and discussion
Simulation with a p53-Mdm2 Network
Truth table of a PBN for the p53-Mdm2 network: x _{ 1 } , x _{ 2 } are the present states of p53 and Mdm2
x _{1} x _{2} | ${\mathit{f}}_{\mathbf{1}}^{\left(\mathbf{1}\right)}$ | ${\mathit{f}}_{\mathbf{2}}^{\left(\mathbf{1}\right)}$ | ${\mathit{f}}_{\mathbf{3}}^{\left(\mathbf{1}\right)}$ | ${\mathit{f}}_{\mathbf{4}}^{\left(\mathbf{1}\right)}$ | ${\mathit{f}}_{\mathbf{1}}^{\left(\mathbf{2}\right)}$ | ${\mathit{f}}_{\mathbf{2}}^{\left(\mathbf{2}\right)}$ | ${\mathit{f}}_{\mathbf{3}}^{\left(\mathbf{2}\right)}$ | ${\mathit{f}}_{\mathbf{4}}^{\left(\mathbf{2}\right)}$ |
---|---|---|---|---|---|---|---|---|
00 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 1 |
01 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |
10 | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 |
11 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 0 |
$\phantom{\rule{0.12em}{0ex}}{\mathit{c}}_{\mathit{j}}^{\left(\mathit{i}\right)}$ | 0.5 | 0.4 | 0.09 | 0.01 | 0.5 | 0.4 | 0.09 | 0.01 |
An internal entry indicates whether a logical 1 or 0 would result at the next state of a gene from a selected Boolean function. $\phantom{\rule{0.12em}{0ex}}{\mathit{c}}_{\mathit{j}}^{\left(\mathit{i}\right)}$ indicates the transition probability by a function ${\mathit{f}}_{\mathit{j}}^{\left(\mathit{i}\right)}$[24].
The network function and selection probability for each context in the p53-Mdm2 network
Context | S2S3, S0S1 | Combination | Selection probability | Context | S2S3, S0S1 | Combination | Selection probability |
---|---|---|---|---|---|---|---|
1 | 00,00 | ${\mathit{f}}_{1}^{\left(2\right)}{\mathit{f}}_{1}^{\left(1\right)}$ | 0.25 | 9 | 10,00 | ${\mathit{f}}_{3}^{\left(2\right)}{\mathit{f}}_{1}^{\left(1\right)}$ | 0.045 |
2 | 00,01 | ${\mathit{f}}_{1}^{\left(2\right)}{\mathit{f}}_{2}^{\left(1\right)}$ | 0.2 | 10 | 10,01 | ${\mathit{f}}_{3}^{\left(2\right)}{\mathit{f}}_{2}^{\left(1\right)}$ | 0.036 |
3 | 00,10 | ${\mathit{f}}_{1}^{\left(2\right)}{\mathit{f}}_{3}^{\left(1\right)}$ | 0.045 | 11 | 10,10 | ${\mathit{f}}_{3}^{\left(2\right)}{\mathit{f}}_{3}^{\left(1\right)}$ | 0.0081 |
4 | 00,11 | ${\mathit{f}}_{1}^{\left(2\right)}{\mathit{f}}_{4}^{\left(1\right)}$ | 0.005 | 12 | 10,11 | ${\mathit{f}}_{3}^{\left(2\right)}{\mathit{f}}_{4}^{\left(1\right)}$ | 0.0009 |
5 | 01,00 | ${\mathit{f}}_{2}^{\left(2\right)}{\mathit{f}}_{1}^{\left(1\right)}$ | 0.2 | 13 | 11,00 | ${\mathit{f}}_{4}^{\left(2\right)}{\mathit{f}}_{1}^{\left(1\right)}$ | 0.005 |
6 | 01,01 | ${\mathit{f}}_{2}^{\left(2\right)}{\mathit{f}}_{2}^{\left(1\right)}$ | 0.16 | 14 | 11,01 | ${\mathit{f}}_{4}^{\left(2\right)}{\mathit{f}}_{2}^{\left(1\right)}$ | 0.004 |
7 | 01,10 | ${\mathit{f}}_{2}^{\left(2\right)}{\mathit{f}}_{3}^{\left(1\right)}$ | 0.036 | 15 | 11,10 | ${\mathit{f}}_{4}^{\left(2\right)}{\mathit{f}}_{3}^{\left(1\right)}$ | 0.0009 |
8 | 01,11 | ${\mathit{f}}_{2}^{\left(2\right)}{\mathit{f}}_{4}^{\left(1\right)}$ | 0.004 | 16 | 11,11 | ${\mathit{f}}_{4}^{\left(2\right)}{\mathit{f}}_{4}^{\left(1\right)}$ | 0.0001 |
The CSSBN is then used to obtain the state transition matrix (STM) for the p53-Mdm2 network. See Additional file 3: The Matlab program that describes the structure of the CSSBN in Figure 6 and computes its STM for the p53-Mdm2 network with perturbation.
Differences in the state transition matrices (STMs) obtained using the context-sensitive stochastic Boolean network (CSSBN) with perturbation, compared to the results by using the analytical context-sensitive probabilistic Boolean network (CSPBN) approach in[23]
L(bits) | q = 1p = 0 | q = 0.99p = 0 | q = 0.5p = 0 | q = 0.8p = 0.01 | |
---|---|---|---|---|---|
||Δ A||_{1} | 1k | 0.1820 | 0.2067 | 0.2320 | 0.3469 |
10k | 0.0754 | 0.0754 | 0.0605 | 0.0890 | |
100k | 0.0291 | 0.0239 | 0.0306 | 0.0210 | |
||Δ A||_{2} | 1k | 0.0642 | 0.0754 | 0.0706 | 0.1012 |
10k | 0.0258 | 0.0259 | 0.0207 | 0.0297 | |
100k | 0.0101 | 0.0091 | 0.0097 | 0.0080 | |
||Δ A||_{∞} | 1k | 0.0382 | 0.0413 | 0.0530 | 0.0888 |
10k | 0.0144 | 0.0145 | 0.0164 | 0.0249 | |
100k | 0.0050 | 0.0055 | 0.0064 | 0.0074 | |
Average time (s) | 1k | 0.0601 | 0.0547 | 0.0534 | 0.0623 |
10k | 0.3172 | 0.3180 | 0.3863 | 0.3186 | |
100k | 2.8781 | 2.8531 | 2.8467 | 2.9226 |
As revealed in Table 4, the difference in the STMs computed using the CSSBN and an analytical CSPBN approach is significantly reduced by increasing the sequence length L. However, the inaccuracies, due to the inherent stochastic fluctuations in stochastic computation, are generally small and thus negligible. Hence, the proposed CSSBN model can be used to accurately and efficiently compute the STM of a CSPBN.
p = 0.01 | p = 0.1 | p = 0.3 | |||
---|---|---|---|---|---|
q = 0.9 | q = 0.8 | q = 0.9 | |||
||Δ SSD _{ AP - AC }||_{1} | 0.0228 | 0.0214 | 0.0041 | ||
||Δ SSD _{ CSSBN - AC }||_{1} | L (bits) | 10k | 0.0094 | 0.0099 | 0.0073 |
50k | 0.0061 | 0.0060 | 0.0057 | ||
100k | 0.0047 | 0.0056 | 0.0040 | ||
||Δ SSD _{ AP - AC }||_{2} | 0.0118 | 0.0118 | 0.0025 | ||
||Δ SSD _{ CSSBN - AC }||_{2} | L (bits) | 10k | 0.0055 | 0.0057 | 0.0041 |
50k | 0.0036 | 0.0036 | 0.0030 | ||
100k | 0.0028 | 0.0031 | 0.0025 | ||
||Δ SSD _{ AP - AC }||_{∞} | 0.0074 | 0.0081 | 0.0020 | ||
||Δ SSD _{ CSSBN - AC }||_{∞} | L (bits) | 10k | 0.0047 | 0.0045 | 0.0027 |
50k | 0.0030 | 0.0029 | 0.0021 | ||
100k | 0.0020 | 0.0022 | 0.0017 |
It has been shown that a PBN with perturbation evolves as an ergodic Markov chain [12]. For a larger perturbation probability p, there is an increased randomness in the network activities, thus the steady states of the network are more evenly distributed [23]. This can be seen in the simulation results in Figure 9. Figures 7 and 9 further show that context switching has little impact on the SSD of the network with perturbation. This is due to the fact that the context switching activity only affects the selection probability of a context, but not the predictor functions. Also shown in Figure 9 is that an oscillation of the states of p53 and Mdm2 exists, indicated by the higher probabilities of the network states (or GAPs) 01 and 10 compared with the other states. However, this two-state oscillation is less evident when the perturbation rate is significantly higher. On the other hand, a modest perturbation to the network only has a minor effect on the oscillatory behavior, confirming the stability and robustness of the p53-Mdm2 network in a dynamic environment [47].
Experiments on a glioma network
F _{ 1 } | F _{ 2 } | F _{ 3 } | F _{ 4 } | F _{ 5 } | F _{ 6 } | F _{ 7 } |
0.8560 | 0.2768 | 0.6759 | 1.0000 | 1.0000 | 0.0263 | 1.0000 |
0.1440 | 0.7232 | 0.3241 | 0.4983 | |||
0.4754 | ||||||
F _{ 8 } | F _{ 9 } | F _{ 10 } | F _{ 11 } | F _{ 12 } | F _{ 13 } | F _{ 14 } |
0.0857 | 1.0000 | 1.0000 | 0.8508 | 1.0000 | 0.8697 | 0.6004 |
0.9143 | 0.1492 | 0.1303 | 0.3996 |
In Table 6, each column for F _{ i }, i ∈ {1, 2, …, 14}, contains the selection probabilities for the Boolean functions of gene i, with the value in the j th row as the probability for ${\mathit{f}}_{\mathit{j}}^{\left(\mathit{i}\right)}$, thus the probabilities sum to 1 in each column. Table 6 also shows that six genes have only a single predictor function; only the state of gene 6 is determined by three functions, while the states of the other genes are determined by two predictor functions.
Steady state analysis
Norms of the differences in the computed steady state distributions (SSDs) and average run time for the glioma network
p = 0.01q = 0.9 | p = 0.1q = 0.9 | ||||||
---|---|---|---|---|---|---|---|
L (bits) | 50k | 200k | 1M | 50k | 200k | 1M | |
||Δ SSD _{ AP - CSSBN }||_{1} | 0.2314 | 0.1323 | 0.0861 | 0.4251 | 0.2104 | 0.0939 | |
||Δ SSD _{ AP - CSSBN }||_{2} | 0.0055 | 0.0032 | 0.0023 | 0.0045 | 0.0023 | 0.0010 | |
||Δ SSD _{ AP - CSSBN }||_{∞} | 0.0013 | 8.01 × 10^{- 4} | 4.13 × 10^{- 4} | 3.05 × 10^{- 4} | 2.16 × 10^{- 4} | 7.85 × 10^{- 5} | |
Average time (s) | CSSBN | 10.106 | 311.55 | 1377.5 | 10.069 | 315.65 | 1380.7 |
Approximate | 21421 | 21417 |
As shown in Table 7, the CSSBN time-frame expansion technique efficiently evaluates the SSD of the glioma network and produces results comparable to those obtained by the approximate analysis [20]. Evaluation accuracy is further improved by using longer stochastic sequences with yet a shorter runtime compared to the approximate method. Since it is difficult, if not impossible, to compute the STM or SSD of a large GRN by using an accurate or approximate analysis, the CSSBN time-frame expansion approach provides an alternative means to accurately and efficiently estimate the SSD of a large network with a tunable accuracy by using stochastic sequences of different lengths.
Intervention analysis
In an n- gene network, if gene X _{ i } is the target gene, the states of all genes can be divided into a set of desirable states, D, and a set of undesirable states, U, by the expression level of the target gene [40]. As can be seen in Figure 10, gene X _{14}, the (BCL2A1); BFL1 protein; GRS protein is one of the most influential genes that interacts closely with others, thus it is chosen as the target gene. The desirable and undesirable states are then separated by the expression level of X _{14}. Assume that the inactivation of X _{14} is preferred; the cumulative probabilities of the desirable and undesirable states are given by $\sum {{{}_{\mathit{x}}}_{{}_{14}}}_{=0}{\mathit{\pi}}_{\mathit{x}}$ and $\sum {{{}_{\mathit{x}}}_{{}_{14}}}_{=1}{\mathit{\pi}}_{\mathsf{y}}$, respectively, where π _{ x } and π _{ y } are elements in the desirable and undesirable SSDs respectively.
Cumulative distributions of the desirable states with a different gene selected as the control gene for the simplified 12-gene context-sensitive glioma network, with perturbation rate p = 0.001 and context switching probability q = 0.99
Gene | X _{1} | X _{2} | X _{3} | X _{4} | X _{6} | X _{8} | X _{9} | X _{10} | X _{11} | X _{12} | X _{13} |
---|---|---|---|---|---|---|---|---|---|---|---|
$\sum _{{\mathrm{x}}_{14}=0}}{\mathrm{\pi}}_{\mathrm{x}$ | 70.91% | 62.93% | 52.61% | 45.90% | 48.98% | 54.78% | 53.76% | 61.61% | 57.23% | 56.99% | 61.78% |
When a different gene is selected as the control gene, the effect of intervention varies with respect to improving the probabilities of desirable states in the SSD. For the glioma network, as revealed in Table 8, gene X_{1}, Tie-2, is the most effective for maximizing the percentage of the desirable states. The cumulative distribution of the desirable states is increased from 45.97% to 70.91% by an intervention using Tie-2 as the control gene. Also revealed in the table is that an intervention via gene X _{4}, (HSP40); (HDJ1; DNAJ1), has nearly no impact on the distribution of the desirable states, as the percentage of the desirable states has not been changed much by the intervention (i.e., 45.90% vs. 45.97%). A modest improvement in the desirable state distribution is obtained by an intervention with another gene as the control gene. For any target gene, this process can be applied to find the most significant gene that maximizes the probabilities of desirable states.
Conclusion
Context-sensitive stochastic Boolean networks (CSSBNs) are proposed as an efficient approach to modeling the effects of gene perturbation and intervention in gene regulatory networks (GRNs). In a CSSBN, the state transition matrix can be accurately and efficiently computed with a complexity of O(nLk 2^{ n }), where n is the number of genes in a context-sensitive probabilistic Boolean network (CSPBN), k is the number of contexts and L is a factor determined by the stochastic sequence length. This result is an improvement compared to the previous result of O(nk ^{2}2^{2n }) for an accurate analysis. The use of non-Bernoulli stochastic sequences further increases computational efficiency and allows for a tunable tradeoff between accuracy and efficiency. A steady state analysis using a time-frame expansion technique has shown a significant speedup and produced more accurate results than an approximate analysis in the computation of the steady state distribution.
CSSBNs are constructed for the analysis of gene perturbation in a p53-Mdm2 network and gene intervention in a glioma network. Biologically meaningful insights are gained into the oscillatory dynamics of the p53-Mdm2 network in a context-switching environment with random gene perturbation. It has been shown that the steady state distribution (SSD) changes drastically with the increase of the perturbation rate, whereas the effect of context switching is rather limited for a given perturbation rate. Therefore, random gene perturbation may have a greater effect on the final distribution of the steady state compared to context switching activities. By predicting the SSD, the CSSBN approach can further evaluate the effectiveness of external gene intervention. A case study of the glioma network shows that the CSSBNs are useful in modeling the effects of gene perturbation and intervention in a complex context-sensitive GRN. This will eventually help drug discovery for an effective drug intervention therapy.
Declarations
Acknowledgements
This work was supported in part by the Natural Sciences and Engineering Research Council (NSERC) of Canada in a Discovery Grant and a Startup fund of the University of Alberta.
Authors’ Affiliations
References
- Elowitz MB, Levine AJ, Siggia ED, Swain PS: Stochastic gene expression in a single cell. Sci. 2002, 297: 1183-1186. 10.1126/science.1070919.View Article
- Pedraza JM, van Oudenaarden A: Noise propagation in genetic networks. Sci. 2005, 307: 1965-1969. 10.1126/science.1109090.View Article
- Dunlop M, Cox RIII, Levine J, Murray R, Elowitz M: Regulatory activity revealed by dynamic correlations in gene expression noise. Nat Genet. 2008, 40: 1493-1498. 10.1038/ng.281.PubMed CentralView ArticlePubMed
- Kauffman SA: Metabolic stability and epigenesis in randomly constructed genetic nets. Theor Biol. 1969, 22: 437-467. 10.1016/0022-5193(69)90015-0.View Article
- Klipp E: Systems Biology in Practice: Concepts, Implementation and Application. 2005, Weinheim, Germany: Wiley-VCHView Article
- Qian L, Wang H, Dougherty E: Inference of noisy nonlinear differential equation models for gene regulatory networks using genetic programming and Kalman filtering. IEEE Trans Signal Process. 2008, 56 (7): 3327-3339.View Article
- Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J Comput Phys. 1976, 22: 403-10.1016/0021-9991(76)90041-3.View Article
- Gillespie DT: Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977, 81: 2340-2361. 10.1021/j100540a008.View Article
- Davidich M, Bornholt S: Boolean network model predicts cell cycle sequence of Fission yeast. PLoS One. 2008, 3 (2): e1672-e2008. 10.1371/journal.pone.0001672.PubMed CentralView ArticlePubMed
- Giacomantonio CE, Goodhill GJ: A Boolean model of the gene regulatory network underlying Mammalian cortical area development. PLoS Comput Biol. 2010, 6 (9): e1000936-10.1371/journal.pcbi.1000936.PubMed CentralView ArticlePubMed
- Huang S: Gene expression profiling, genetic networks, and cellular states: an integrating concept for tumorigenesis and drug discovery. J Mol Med. 1999, 77: 469-480. 10.1007/s001099900023.View ArticlePubMed
- Shmulevich I, Dougherty ER, Zhang W: From Boolean to probabilistic Boolean networks as models of genetic regulatory networks. Proc IEEE. 2002, 90: 1778-1792. 10.1109/JPROC.2002.804686.View Article
- Dougherty ERSI, Zhang W: Gene perturbation and intervention in probabilistic Boolean networks. Bioinform. 2002, 18 (10): 1319-1331. 10.1093/bioinformatics/18.10.1319.View Article
- Shmulevich I, Dougherty E: Probabilistic Boolean Networks: The Modeling and Control of Gene Regulatory Networks. 2009, Philadelphia: SIAM
- Pal R, Datta A, Bittner ML, Dougherty ER: Intervention in context-sensitive probabilistic Boolean networks. Bioinform. 2005, 21: 1211-1218. 10.1093/bioinformatics/bti131.View Article
- Shmulevich I, Gluhovsky I, Hashimoto RF, Dougherty ER, Zhang W: Steady-state analysis of genetic regulatory networks modeled by probabilistic Boolean networks. Comp Funct Genomics. 2003, 4: 601-608. 10.1002/cfg.342.PubMed CentralView ArticlePubMed
- Rosenthal JS: Minorization conditions and convergence rates for Markov chain Monte Carlo. J Am Stat Assoc. 1995, 90: 558-566. 10.1080/01621459.1995.10476548.View Article
- Zhang S, Ching W, Ng M, Akutsu T: Simulation study in probabilistic Boolean network models for genetic regulatory networks. Int J Data Min Bioinform. 2007, 1 (3): 217-240. 10.1504/IJDMB.2007.011610.View ArticlePubMed
- Brun M, Dougherty ER, Shmulevich I: Steady-state probabilities for attractors in probabilistic Boolean networks. Signal Process. 2005, 85: 1993-2013. 10.1016/j.sigpro.2005.02.016.View Article
- Ching W, Zhang S, Ng M, Akutsu T: An approximation method for solving the steady-state probability distribution of probabilistic Boolean networks. Bioinform. 2007, 23: 1511-1518. 10.1093/bioinformatics/btm142.View Article
- Datta A, Choudhary A, Bittner ML, Dougherty ER: External control in markovian genetic regulatory networks. Mach Learn. 2003, 52: 169-191. 10.1023/A:1023909812213.View Article
- Datta A, Choudhary A, Bittner ML, Dougherty ER: External control in markovian genetic reegulatory networks: the imperfect information case. Bioinform. 2004, 20: 924-993. 10.1093/bioinformatics/bth008.View Article
- Pal R: Context-sensitive probabilistic Boolean networks: steady-state properties, reduction, and steady-state approximation. IEEE Trans Signal Process. 2010, 58 (2): 879-890.View Article
- Liang J, Han J: Stochastic Boolean networks: an efficient approach to modeling gene regulatory networks. BMC Syst Biol. 2012, 6: 113-10.1186/1752-0509-6-113.PubMed CentralView ArticlePubMed
- Lahav G, Rosenfeld N, Sigal A, Geva-Zatorsky N, Levine AJ, Elowitz MB, Alon U: Dynamics of the p53-Mdm2 feedback loop in individual cells. Nat Genet. 2004, 36: 147-150. 10.1038/ng1293.View ArticlePubMed
- Martin S, Zhang Z, Martino A, Faulon J-L: Boolean dynamics of genetic regulatory networks infer red from microarray time series data. Bioinform. 2007, 23 (7): 866-874. 10.1093/bioinformatics/btm021.View Article
- Zhu P, Han J: Stochastic multiple-valued gene networks. IEEE Transac Biomed Circ Syst. 2014, 8 (1): 42-53.View Article
- Zhu P, Han J: Asynchronous stochastic Boolean networks as gene network models. J Comput Biol. 2014, in press
- Faryabi B, Vahedi G, Chamberland JF, Datta A, Dougherty ER: Intervention in context-sensitive probabilistic Boolean networks revisited. EURASIP J Bioinform Syst Biol. 2009, 2009: 5-View Article
- McAdams HH, Shapiro L: Circuit simulation of genetic networks. Sci. 1995, 269 (5224): 650-10.1126/science.7624793.View Article
- Abdi A, Tahoori MB, Emamian ES: Fault diagnosis engineering of digital circuits can identify vulnerable molecules in complex cellular pathways. Sci Signal. 2008, 1 (42): ra10-View ArticlePubMed
- Adar R, Benenson Y, Linshiz G, Rosner A, Tishby N, Shapiro E: Stochastic computing with biomolecular automata. PNAS. 2004, 101 (27): 9960-9965. 10.1073/pnas.0400731101.PubMed CentralView ArticlePubMed
- Benenson Y, Gil B, Ben-Dor U, Adar R, Shapiro E: An autonomous molecular computer for logical control of gene expression. Nat. 2004, 429: 423-429. 10.1038/nature02551.View Article
- Gaines BR: Stochastic computing systems. Adv Inf Syst Sci. 1969, 2: 37-172.View Article
- Han J, Chen H, Liang J, Zhu P, Yang Z, Lombardi F: A Stochastic Computational Approach for Accurate and Efficient Reliability Evaluation. IEEE Transactions on Computers. 2013, in press. Advance access in IEEE Xplore
- Pal R, Datta A, Bittner ML, Dougherty ER: Optimal infinite horizon control for probabilistic Boolean networks. IEEE Trans Signal Process. 2006, 54: 2375-2387.View Article
- Faryabi B, Datta A, Dougherty ER: On approximate stochastic control in genetic regulatory networks. IET Syst Biol. 2007, 1 (6): 361-368. 10.1049/iet-syb:20070015.View ArticlePubMed
- Layek R, Datta A, Pal R, Dougherty ER: Adaptive intervention in probabilistic Boolean networks. Bioinform. 2009, 25 (16): 2042-2048. 10.1093/bioinformatics/btp349. Doi:View Article
- Qian X, Dougherty ER: Effect of function perturbation on the steady-state distribution of genetic regulatory networks: optimal structural intervention. IEEE Trans Signal Process. 2008, 56: 4966-4976.View Article
- Qian X, Ivanov I, Ghaffari N, Dougherty ER: Intervention in gene regulatory networks via greedy control policies based on long-run behavior. BMC Syst Biol. 2009, 3: 16-10.1186/1752-0509-3-16.View Article
- Weinberg RA: The Biology of Cancer. 2006, New York: Garland Science, 1
- Vogelstein B, Lane D, Levine AJ: Surfing the p53 network. Nat. 2000, 408: 307-310. 10.1038/35042675.View Article
- Batchelor E, Loewer A, Lahav G: The ups and downs of p53: understanding protein dynamics in single cells. Nat Rev Cancer. 2009, 9 (5): 371-377. 10.1038/nrc2604.PubMed CentralView ArticlePubMed
- Ciliberto A, Novak B, Tyson JJ: Steady states and oscillations in the p53-Mdm2 network. Cell Cycle. 2005, 4: 486-493.
- Abou-Jaoude W, Ouattara D, Kaufman M: From structure to dynamics: frequency tuning in the p53-mdm2 network: I: logical approach. J Theor Biol. 2009, 258 (4): 561-577. 10.1016/j.jtbi.2009.02.005.View ArticlePubMed
- Murrugarra D, Veliz-Cuba A, Aguilar B, Arat S, Laubenbacher R: Modeling stochasticity and variability in gene regulatory networks. EURASIP J Bioinform Syst Biol. 2012, 1: 5-View Article
- Fuller GN, Rhee CH, Hess KR, Caskey LS, Wang R, Bruner JM, Yung WK, Zhang W: Reactivation of insulin-like growth factor binding protein 2 expression in glioblastoma multiforme: a revelation by parallel gene expression profiling. Cancer Res. 1999, 59: 4228-4232.PubMed
- Ge H, Qian M: Boolean network approach to negative feedback loops of the p53 pathways: synchronized dynamics and stochastic limit cycles. J Comput Biol. 2009, 16: 119-132. 10.1089/cmb.2007.0181. doi:10.1089/cmb.2007.0181View ArticlePubMed
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.