2.1. Definition

A context-sensitive probabilistic Boolean network consists of a sequence V = {x_i}_{i = 1}ⁿ of n nodes, where x_i ∈ {0,…, d − 1}, and a sequence {f_l}_{l = 1}^k of vector-valued functions defining constituent networks. In the framework of gene regulation, each element x_i represents the expression value of a gene. It is common to abuse terminology by referring to x_i as the ith gene. Each vector-valued function f_l = (f_l1,…, f_ln) represents a constituent network of the context-sensitive PBN. The function f_li : {0,…, d − 1}ⁿ → {0,…, d − 1} is the predictor of gene i, whenever context l is selected. The number of quantization levels for gene expressions is denoted by d. At each updating epoch, a random variable determines whether the constituent network is switched or not. The switching probability q is a system parameter. If the context remains unchanged, then the context-sensitive PBN behaves like a fixed Boolean network where the values of all the genes are updated synchronously according to the current constituent network. On the other hand, if a switch occurs, then a constituent network is randomly selected from {f_l}_{l = 1}^k according to the selection probability distribution {r_l}_{l = 1}^k. Once the predictor function f_l is determined, the values of the genes are updated using the new constituent network, that is, according to the rules defined by f_l.

Two quantization levels have thus far been used in practice. If d = 2 (binary), then the constituent networks are Boolean networks with 0 or 1 meaning OFF or ON, respectively [¹⁰]. The case where d = 3 (ternary) arises when we consider individual genes to be downregulated (0), upregulated (2), or invariant (1). This situation commonly occurs with cDNA microarrays, where a ratio is taken between the expression values on the test channel (red) and the base channel (green) [¹⁶]. In this paper, we will develop the methodology for d = 2, so that gene values are either 0 or 1. The methodology can be extended to other finite quantization levels, albeit, at the expense of tedious mathematical expressions. All the binary operations in this section would need to be replaced by case statements, and the perturbation process should be articulated on a case by case basis.

We focus on context-sensitive PBNs with perturbations, meaning that each gene may change its value with small probability p at each epoch. If γ_i(t) is a Bernoulli random variable with parameter p and the random vector γ at instant t is defined as γ(t) = (γ₁(t), γ₂(t),…, γ_n(t)), then the value of gene i is determined at each epoch t by

where operator ⊕ is componentwise addition in modulo two and f_li is the predictor of gene i according to the current context of the network l. Such a perturbation captures the realistic situation where the activity of a gene is subject to random alteration.

The gene-activity profile (or GAP) is an n-digit binary vector x(t) = (x₁(t),…, x_n(t)) giving the expression values of the genes at time t, where x_i(t) ∈ {0, 1}. We denote the set of all possible GAPs by 𝒳. The dynamic behavior of a context-sensitive PBN can be modeled by a Markov chain whose states are ordered pairs consisting of a constituent network κ and a GAP x. The evolution of the context-sensitive PBN can therefore be represented using a stationary discrete-time equation

where state z(t) is an element of the state space 𝒵 = {(κ, x) : κ ∈ {1,…, k}, x ∈ 𝒳}. The disturbance w(t) is the manifestation of uncertainties in the biological system, due either to context switching or a change in the GAP resulting from a random gene perturbation. It is assumed that both the gene perturbation distribution and the network switching distribution are independent and identically distributed over time. The switching frequency of the context differentiates the instantaneously random PBN from the context-sensitive PBN. If the contexts change at every instant, that is, q = 1, then the PBN is instantaneously random.

We note that a bijection can be drawn between the gene-activity profile x(t) or the states z(t) and their decimal representations x(t) and z(t) based on their binary expansion. The integers x(t) and z(t) take values in X = {0, 1,…, dⁿ − 1} and Z = {0, 1,…, k × dⁿ − 1}, respectively. These decimal representations facilitate the depiction of our numerical results in Section 3.

2.2. Infinite-Horizon Intervention

We can formulate the task of finding the most effective intervention strategy as a sequential decision making problem, when the dynamics of a context-sensitive PBN are expressed according to (2). To this end, we can specify the Markov chain that describes the dynamics of the context-sensitive PBN by defining its transition probability matrix and initial state distribution.

In the presence of an external regulator, we suppose that the context-sensitive PBN has a binary intervention input u_g(t) on the control gene g. The intervention input u_g(t), which takes values in set 𝒞 = {0, 1}, specifies the action on the control gene g. Treatment alters the status of the control gene g, which can be selected from all the genes in the network. If treatment is applied, u_g(t) = 1, then the state of the control gene g is toggled; otherwise the state of the control gene g remains unchanged.

For the case of a single control gene g, the system evolution is represented by a stationary discrete-time equation

where the state z(t) is an element of 𝒵, and similar to the context-sensitive PBN without control, w(t) is the manifestation of uncertainties in the model. The transition probability matrix for the controlled context-sensitive PBN can be defined easily, once the transition probability matrix of the uncontrolled system is known. We derive an expression for this matrix in Section 2.3. Originating from a state z₁, the successor state z₂ is selected randomly within set 𝒵 according to the transition probability P(z₁, z₂, u);for all z₁, z₂ ∈ 𝒵 and all u ∈ 𝒞.

To define the problem of intervention in a context-sensitive PBN, we associate a cost-per-stage c(z₁, z₂, u) to each possible event. In general, the cost-per-stage can depend on the origin state z₁, the successor state z₂, and the control input u. We define the average immediate cost in state z₁, when control u is selected, by

We consider a discounted formulation of the expected total cost. The discounting factor, λ ∈ (0, 1), ensures convergence of the expected total cost over the long run [¹⁷]. In the case of cancer therapy, the discounting factor also emphasizes that obtaining treatment at an earlier stage is favored over later stages.

For initial state i and strategy π_g = {μ_g(·, 0), μ_g(·, 1),…}, where μ_g(·, t) : 𝒵 → 𝒞 denotes the decision rule at epoch t, the infinite expected total discounted cost is defined by

The sequential decision maker must identify an optimal strategy π_g^∗ = {μ_g^∗(·, 0), μ_g^∗(·, 1),…} such that the J_πg(i) is minimized for each state i ∈ 𝒵. Mathematically, an optimal strategy π_g^∗ is a solution of the optimization problem

For the specifics of our formulation, an optimal strategy always exists [¹⁷]. It is given by the fixed-point solution of the Bellman optimality equation

Moreover, this optimal strategy is stationary and independent of the initial state i [¹⁷]. Standard dynamic programming algorithms can be used to find a fixed-point of the Bellman optimality equation. In our model, gene perturbation ensures that all the states communicate with one another. Hence, the Markov decision process associated with any stationary policy is ergodic and has a unique invariant distribution equal to its limiting distribution [¹⁸].

2.3. Transition Probability Matrix of a Context-Sensitive Probabilistic Boolean Network

For a given cost of intervention and cost-per-stage, the solution to (8) depends on the transition probability matrix in (4). The latter can be found by observing that two mutually exclusive events may occur at any epoch: the current context of the network remains the same for two consecutive instants, or the context of the network changes to a new one at the time instant t + 1. Moreover, the context may remain unchanged in two mutually exclusive ways: the binary switching variable is 0, which means that no change is possible, or the binary switching variable is 1, and the current network is picked again through random selection [¹¹, ¹²]. In particular, when the switching variable is 1, a new context is selected independent of the current system state. Thus, the same network can be active twice in a row. This interpretation of switching the context in a PBN is in concert with the original definition of context-sensitive PBNs in [¹²]. Before proceeding, we note that transitioning was defined differently in [¹³, ¹⁴], where it was assumed that, whenever the switching variable is 1, a change of context must occur; the result being that context selection is conditioned on the current context. While this contrast produces a difference in the transition probabilities, it does not change the underlying issue in this paper, that being analyzing the effects of the state-space reduction proposed in [¹³, ¹⁴] on the performance of therapeutic interventions.

Letting z₁ = (κ₁, x₁) and z₂ = (κ₂, x₂) be two states in 𝒵, we derive the transition probability

from z₁ to z₂ in the absence of intervention. Note that we can rewrite expression (9) asUsing the Bayes theorem, we getwhere 1(·) is the indicator function. Furthermore, we haveand when κ₁ ≠ κ₂, we getHere, q is the probability of switching context, and r_κi is the probability of selecting context κ_i during a switch.

A transition from GAP x₁ to GAP x₂ may occur either according to the constituent network at instant t + 1 or through an appropriate number of random perturbations, but not both. Let us define F_l by

Then, we haveand, for κ₁ ≠ κ₂, we obtainwhere D(x₁, x₂) is the Hamming distance between two gene-activity profiles x₁ and x₂.

The first parts of (15) and (16) correspond to the probability of transition from GAP x₁ to GAP x₂ according to the predictor functions defined by the constituent network at time instant t + 1. The remaining terms account for transition between GAPs that are due to random gene perturbation.

By replacing the terms of expression (11) by their equivalents from (12), (13), (15), and (16), it can be shown that the probability of transition from any state z₁ = (κ₁, x₁) to z₂ = (κ₂, x₂) is given by

The elements of the transition probability matrix of the controlled context-sensitive PBN given by (4) can then be expressed through (17). The value of state after intervention z₁ = (κ₁, x₁) can be determined according to the status of the control signal and the value of the state prior to the intervention z^1=(κ^1,x^1). Here, κ₁ is equated to κ^1, and the value of the GAP is updated according to the value of the control signal in the devised strategy according to

All the 2ⁿ elements of vector e_g are zeros except the element at the gth position, which is set to one.

2.4. Approximate Transition Probability Matrix of a Context-Sensitive Probabilistic Boolean Network

Following the reduction method proposed in [¹³, ¹⁴], we derive an expression for the approximate transition probability matrix in which the context is removed from the state space of the system. We base our derivation on the stochastic matrix defined by (17). The approximate stochastic model describes the dynamics of the system solely based on the GAPs, and its state space takes values from the set 𝒳. The probability of transition from GAP x₁ to GAP x₂ in two consecutive epochs is derived from the weighted sum of the actual transition probabilities with respect to the selection probabilities of the contexts. If we denote the probability of transition between two GAPs by

then under the reduction assumptions, we defineMoreover, we can expand this expression aswhich in turn can be presented aswhereThe above expression for Λ can be further simplified asThus, we haveThe last expression for Λ can be further reduced toby setting ∑_κ1r_κ1 = 1.

Equations (22) and (26) express the approximate transition probability matrix associated with the reduced model. Although we have started from a different expression for the transition probability matrix of a context-sensitive PBN owing to a different interpretation of switching contexts, our final expression for the approximate transition probability matrix is similar to the one in [¹³, ¹⁴]. Averaging over the various contexts in (20) reduces the transition probability distributions associated with a context-sensitive PBN to transition probability distributions arising from the corresponding instantaneously random PBN, the fact that is overlooked in [¹³, ¹⁴]. The transition probability matrix of the corresponding instantaneously random PBN with the similar parameters can be obtained from expression (17), when the context is allowed to switch at each epoch by setting q = 1. Hence, the optimal and approximate intervention strategies perform similarly whenever the switching probability approaches to value one. This observation is supported by our numerical studies.

3.1. Synthetic Networks

In our numerical studies, we postulate the following cost-per-stage:

where c is the cost of control, and 𝒰, 𝒟 are the sets of undesirable and desirable states, respectively. We set c = 1 to make the application of intervention more plausible compared to visiting undesirable states. The target gene is chosen to be the most significant gene in the GAP. We assume that the upregulation of the target gene is undesirable. Consequently, the state space is partitioned into desirable states, 𝒟 = {(κ, x) : κ ∈ {1,…, k}, x ∈ {0,…, 2ⁿ⁻¹ − 1}}, and undesirable states, 𝒰 = {(κ, x) : κ ∈ {1,…, k}, x ∈ {2ⁿ⁻¹,…, 2ⁿ − 1}}, where n is the number of genes. We use the natural decimal bijection of the GAP x to facilitate the presentation of our results. We set the number of genes to be five. The study of networks with larger numbers of genes would be computationally prohibitive due to the complexity of the corresponding dynamic program. The cost values have been chosen in accordance with an earlier study [¹⁴]. Since our objective is to downregulate the target gene, a higher cost is assigned to destination states having an upregulated target gene. Moreover, for a given status of the target gene, a higher cost is assigned when the control is applied, versus when it is not. In practice, the cost values have to mathematically capture the benefits and costs of intervention and the relative preference of states. They must be set with the help of physicians in accordance with their clinical judgement. Although this is not feasible within the realm of current medical practice, we do believe that such an approach will become feasible when engineering approaches are integrated into translational medicine.

We generate synthetic context-sensitive PBNs in the following manner. Each context-sensitive PBN consists of two contexts. Each constituent network is randomly generated with bias equal to 0.5. The bias is the probability that a randomly generated Boolean function takes on a value of one. To complete the specification of a context-sensitive PBN, we need to specify the selection and switching probability distributions, along with its constituent networks. We consider the selection and switching probabilities as parameters during our numerical study. In the first set of experiments, we assume that the constituent networks are selected with equal probabilities. Then, we vary the value of the switching probability. In the second set of simulations, the switching probability is fixed but the selection probability varies. We generate one thousand random synthetic context-sensitive PBNs for each scenario. Our objective is to utilize the statistics generated by these synthetic networks to evaluate the effects of removing the context from the state space of a context-sensitive PBN when designing an intervention strategy.

For each context-sensitive PBN, the exact and approximate transition probability matrices are computed according to (17) and (22), respectively. Thereafter, we solve the optimal intervention problems for the original model and its reduced approximation and find the corresponding optimal strategies.

The devised strategy, μ_g^∗ : 𝒵 → 𝒞, for the exact transition probability matrix specifies the action that should be taken at each time step. The second policy is based on the reduced stochastic matrix and only takes the GAP as its input. Since the performance of the approximate strategy must be evaluated with respect to the dynamics specified by the original model, we need to extend the approximate strategy to elements of 𝒵. This is achieved by simply disregarding the context element of state z(t) and determining the action based on its GAP element. We denote the resulting intervention strategy obtained through state collapse by μ^g:𝒵 → 𝒞.

In the first set of experiments, we determine the effect of the switching probability q on overall performance. To this end, we generate one thousand context-sensitive PBNs for each value of the switching probability. The selection probability is assumed to have uniform distribution for all the generated context-sensitive PBNs. We set the perturbation probability to 0.01 for all simulations.

For each context-sensitive PBN generated per the above method, we select a random control gene. Then, we use dynamic programming and derive an optimal intervention strategy μ_g^∗ based on the exact transition probability matrix (17). Similarly, an optimal strategy based on the approximate transition probability matrix (22) is derived for the same control gene and is extended to the approximate strategy μ^g. For a context-sensitive PBN, we estimate the average total discounted cost induced by the given optimal strategy μ_g^∗. To this end, we generate synthetic time-course data for thousand time steps from the transition probability matrix for the context-sensitive PBN, while intervening based on optimal strategy μ_g^∗. We estimate the total cost by accumulating the discounted cost of each state given the action at that state. This procedure is repeated ten thousand times for random initial states, and the average of the induced total discounted costs is computed. Following a similar procedure, the approximate strategy μ^g is applied to the system, and the average total discounted cost is computed. Finally, we compute the average total discounted cost for time-course data when no intervention is applied. From here on, we omit the subscript g from the notation of a strategy μ_g to simplify our notations. Since the control gene is selected randomly, this will not affect the following discussions.

The effectiveness of an intervention strategy can be evaluated by computing the difference between its induced cost and the cost accumulated in the absence of intervention. For each set of constituent networks and a given switching probability, we compute the following functions: J¯μ∗, J¯μ^, and J¯. These are the average total discounted cost for a given context-sensitive PBN induced by applying optimal strategy μ^∗, approximate strategy μ^, and no intervention, respectively. The preceding procedure is repeated for one thousand random context-sensitive PBNs, thereby yielding one thousand values for each statistic. We compare the effects of these strategies by computing averages denoted by E[J¯μ∗], E[J¯μ^], and E[J¯].

We consider the percentage of reduction in the average total discounted cost as a performance metric. The normalized gain obtained by each intervention strategy is taken as the immediate consequence of the intervention formulation. This metric is defined as the difference between the average discounted cost before and after intervention, normalized by the cost before intervention. The normalized gain corresponding to the optimal strategy μ^∗ is

and the normalized gain corresponding to the strategy derived from the approximate method μ^ isFigure 1 depicts the results of the first experiment, where q is the parameter of interest. As q increases to one, the difference between normalized gains ΔJ^A and ΔJ^E decreases. The approximating method yields close to optimal performance when the switching probability is large, which is outside the range of typical values used for context-sensitive PBNs. If one cannot obtain context knowledge or the number of contexts results in an unacceptable computational burden, the approximate method provides a strategy for the realistic value q = 0.01, which yields a 30% reduction in performance.

As a byproduct of the intervention formulation, we also consider the effect of an intervention strategy μ on the amount of change in the steady-state probability of undesirable states before and after the intervention. For each set of constituent networks and for a given switching probability, we compute ΔP_μ^∗ and ΔPμ^. These are the normalized reduction in the total probability of visiting undesirable states in the long run for a given context-sensitive PBN when strategies μ^∗ and μ^ are applied to original system, respectively. In other words, we define

where π_μ^∗(i) is the probability of being in state i in the long run under optimal strategy μ^∗; πμ^(i) is the probability of being in state i in the long run under approximate strategy μ^; π(i) is the probability of being in state i in the long run when no control is applied.

The preceding procedure is repeated for one thousand random context-sensitive PBNs, thereby yielding one thousand values for each statistic. We compare the effect of the strategies devised by the exact and approximate transition probability matrices via the empirical averages of each sample sequence, denoted by ΔP^E and ΔP^A. Figure 2 shows ΔP^A and ΔP^E as functions of the switching probability. The trends are similar to those observed for the normalized gains.

In practice, treatment options, such as chemotherapy, have detrimental side effects. A large number of interventions can cause collateral damage that reduces a patient's quality of life. Thus, we define the quantity Γ_μ as the expected number of interventions when the strategy μ is applied in the long run to gauge these side effects. In particular, Γ_μ^∗ and Γμ^ are the expected numbers of executed interventions in the long run using the optimal strategy μ^∗ and the approximate strategy μ^, respectively. We define

where π_μ^∗(i) and πμ^(i) have similar definitions as in (30).

The preceding procedure is repeated for one thousand random context-sensitive PBNs. We compare the expected number of executed interventions using the difference in empirical averages, denoted by ΔΓ≜Γμ^−Γμ∗. Figure 3 indicates the variation in ΔΓ as a function of switching probability q. According to this figure, for small switching probabilities, the approximate strategy μ^ is likely to cause more detrimental side effects.

We study the effect of selection probability on the performance of the approximate strategy μ^ in a second set of experiments. We follow the same procedure as before, except that we set q = 0.01, and we vary the probability of selecting each constituent network. We consider two constituent networks so that the selection probabilities are a function of r₁, the probability of selecting the first context. From Figure 4, as |r₁ − r₂| gets smaller, the difference between the performance of strategies μ^∗ and μ^ diminishes.

Figure 5 compares the steady-state measures ΔP^E and ΔP^A for the optimal and approximate strategies, respectively. The most interesting observation is that, whereas ΔJ^E − ΔJ^A decreases as |r₁ − r₂| decreases, ΔP^E − ΔP^A increases. These different behaviors are not contradictory, since the intervention strategy is designed to minimize the total cost and the improvement in the steady-state behavior is a side effect of our goal. We observe that both ΔJ^E and ΔP^E are stable across parameters, whereas the metrics ΔJ^A and ΔP^A vary considerably. The context removal approximation affects both ΔJ^A and ΔP^A. That is not the case for the exact transition probability matrix. We suspect that a mathematical analysis of this effect is complicated since it involves interaction between the optimization and the reduction. Finally, ΔΓ is plotted as a function of the selection probability in Figure 6. Here, we observe that ΔΓ increases as |r₁ − r₂| decreases.

3.2. A Melanoma Case Study

In this section, we compare the performance of the optimal and approximate strategies in the context of a gene regulatory network developed from steady-state data. This steady-state data was collected in a profiling study of metastatic melanoma in which high abundance of messenger RNA for the gene WNT5A was found to be highly discriminating between cells with properties typically associated with high metastatic competence versus those with low metastatic competence [¹⁹]. Seven genes were considered in [¹³, ¹⁴]: WNT5A, pirin, S 100 P, RET 1, MART 1, HADHB, and STC 3. We apply the design procedure proposed in [²⁰] to generate a context-sensitive PBN possessing four constituent networks. The method of [²⁰] generates Boolean networks with given attractor structures, and the overall context-sensitive PBN is designed so that the data points, which are assumed to come from the steady-state distribution of the network, are attractors in the resulting network. The regulatory graphs of these constituent networks can be found in [¹⁴]. This approach is reasonable because our interest is in controlling the long-run behavior of the network. The intervention objective for this 7-gene network is to downregulate WNT5A. The gene WNT5A ceasing to be downregulated is strongly predictive of the onset of metastasis. A number of other intervention studies based on the same data have aimed to downregulate WNT5A. This model has been used since the discovery of the relation between WNT5A and metastasis. The binary nature of the up or down regulation suits our binary model. A state is desirable, that is, belongs to 𝒟, if WNT5A = 0, and undesirable, that is, belongs to 𝒰, if WNT5A = 1. As we mentioned earlier, application of intervention requires the designation of desirable and undesirable states, and this depends upon the existence of relevant biological knowledge. The use of WNT5A is one such example where the knowledge of practitioners is incorporated in a theoretical framework. Based on our objective, the cost of control is assumed to be one, and the states are assigned penalties according to the cost-per-stage (27). This is the same cost structure as in [¹⁴]. Since our objective is to downregulate WNT5A, a higher penalty is assigned for states having WNT5A upregulated. Also, for a given WNT5A status, a higher penalty is assigned when the control signal is active versus when it is not.

The optimal and approximate intervention strategies are found for the melanoma-related context-sensitive PBN when different genes in the network (except WNT5A itself) are employed as the control genes. Figure 7 depicts the normalized gains when the optimal and approximate strategies for each control gene are used to intervene in the context-sensitive PBN. To compute the normalized gains, we computed the costs for ten thousand trajectories of length two hundred thousand. As we expected, the optimal strategy outperforms the approximate strategy significantly for all the control genes. Moreover, for the best control gene S100P, the difference between the two strategies is the greatest.

Figure 8 depicts the effects of the optimal and approximate strategies on the normalized reduction in the aggregated long-run probability of visiting undesirable states ΔP^E and ΔP^A, respectively. Here, the strategy based on the S100P outperforms the strategies devised for other control genes. Note that the performance differences are not significant for most of the control genes. In particular, one should not draw any conclusions from the fact that ΔP^E is slightly less than ΔP^A in a couple of cases. The intervention strategy is designed to minimize the total cost, and the improvement in the steady-state behavior is a side effect of our method.

Lastly, Figure 9 shows the difference between the expected number of executed interventions for the optimal strategy and the one derived from the approximate representation of the system. Note that the approximate strategy based on the most effective control gene applies 35% more interventions compared to the optimal one, while its performance is still worse.