Exact solving and sensitivity analysis of stochastic continuous time Boolean models

Background Solutions to stochastic Boolean models are usually estimated by Monte Carlo simulations, but as the state space of these models can be enormous, there is an inherent uncertainty about the accuracy of Monte Carlo estimates and whether simulations have reached all attractors. Moreover, these models have timescale parameters (transition rates) that the probability values of stationary solutions depend on in complex ways, raising the necessity of parameter sensitivity analysis. We address these two issues by an exact calculation method for this class of models. Results We show that the stationary probability values of the attractors of stochastic (asynchronous) continuous time Boolean models can be exactly calculated. The calculation does not require Monte Carlo simulations, instead it uses graph theoretical and matrix calculation methods previously applied in the context of chemical kinetics. In this version of the asynchronous updating framework the states of a logical model define a continuous time Markov chain and for a given initial condition the stationary solution is fully defined by the right and left nullspace of the master equation’s kinetic matrix. We use topological sorting of the state transition graph and the dependencies between the nullspaces and the kinetic matrix to derive the stationary solution without simulations. We apply this calculation to several published Boolean models to analyze the under-explored question of the effect of transition rates on the stationary solutions and show they can be sensitive to parameter changes. The analysis distinguishes processes robust or, alternatively, sensitive to parameter values, providing both methodological and biological insights. Conclusion Up to an intermediate size (the biggest model analyzed is 23 nodes) stochastic Boolean models can be efficiently solved by an exact matrix method, without using Monte Carlo simulations. Sensitivity analysis with respect to the model’s timescale parameters often reveals a small subset of all parameters that primarily determine the stationary probability of attractor states.


Comparison of exact method with MaBoSS simulations
The derivation of the stationary states from the full dynamical equations in Section 'Derivation' of the main text (adopted from [1]) means that the values of the stationary probabilities from stochastic simulations and the exact calculation method must be identical, up to a limit due to noise in Monte Carlo simulations that increases with a lower number of sample trajectories. To test if our implementation of the exact calculation method is correct, we compared the results from MaBoSS simulations to calculations in our ExaStoLog toolbox (https://github.com/mbkoltai/exact-stoch-log-mod) for several different logical models.
The full dynamics of model state probabilities with a small 3-node model is shown in Figure 1 We similarly compared results for the cell cycle model of [2] that has a large cyclic attractor of 270 states, again finding discrepancies not larger than 1%. The MaBoSS files for this simulation are model files/maboss/traynard2016 mammalian cellcycle.bnd, model files/maboss/traynard2016 mammalian cellcycle.cfg.
The probability values of states from MaBoSS simulations of a cyclic attractor were averaged over the last 5000 time steps of the simulations to be comparable to the exact solution. We can again see on Figure  2 that the exact method is identical to stochastic simulations up to 1% deviations. We similarly compared ExaStoLog calculations to MaBoSS simulations for other models with fixed point attractors, again getting identical results.
2 List of models analyzed 2.1 [Traynard 2015]: mammalian cell cycle model (13 nodes) Source file at GitHub repository: model files/mammalian cc.bnet This is a 13-node model of the mammalian cell cycle, published in [2], its influence graph shown on SI Fig. 3. The model has two separate subgraphs. In one subgraph there are 2 fixed point attractors (where Rb b2, p27 b2=1 and CycD, Skp2=0, the value of CycA differentiating the two states), in the other a large cyclic attractor of 270 states. We analyzed if the stationary probabilities of the two fixed point attractors and the states within the cyclic attractor are sensitive to the transition rates. Results of parameter scans are shown on SI Fig. 9.

[Zanudo 2017]: breast cancer model (20 nodes)
Source file at GitHub repository: model files/breast cancer zanudo2017.bnet This model is the 20-node subnetwork of the full model published in [3], shown on Figure 5A of [3] in the main text. The influence graph is shown on SI Fig. 4. The model has 4 input nodes: Alpelisib, Everolimus, PDK1, PIM, and the state transition graph is made up of 36 disconnected subgraphs. Due to the high number of disconnected subgraphs if we populate all of them by uniform initial conditions we have a large number of attractor states, shown on SI Fig. 5, all of which are separate stable states. The calculation time (see Table 1 in main text) depends on the choice of initial conditions and how many of the subgraphs are populated. Calculation for one subgraph is below one second. The model has two phenotypic nodes: Proliferation and Apoptosis. We analyze the dependence of the stationary probabilities of these two variables under different initial conditions in terms of the value of Alpelisib, Everolimus and the PIM node, as shown on Figure 4 of the main text. The initial values of these input nodes can be interpreted as the presence of the two drugs and the mutational status of the PIM family of oncoproteins. There are two input nodes: DNAdamage and ECMicroenv. We analyzed the parameter-dependence of the attractor states, characterized by co-activation of one group of phenotypic nodes under the 4 different initial conditions. The results of the sensitivity analysis are summarized on Figures 5, 6 of the main text as well as SI Fig. 10. We have learned from the sensitivity analysis that it is primarily the transition rates d miRNA, u p53, d p63 73, u EMTreg and d Notch pthw that determine the cell fate decision between apoptosis and EMT/proliferation/invasion. Also, we have recapitulated a known (see [4]) synergistic effect of p53 loss-offunction and Notch pathway gain-of-function mutations in a continuous sense by our two-dimensional scan shown on SI Fig. 11.

[Sahin 2009]: breast cancer ERBB-model (20 nodes)
Source file at GitHub repository: model files/sahin breast cancer refined.bnet This is a 20-node model on the role of ERBB receptor overexpression/mutations and its inhibition by trastuzumab (a monoclonal antibody) in breast cancer from [5]. The influence graph is shown on SI Fig.  7. The model's output node is pRB that represents the phosphorylation of the retinoblastoma protein, that is a proxy for G1/S transition (cell cycle progression) and proliferation. This model's attractor is robust to   3 Parameter sensitivity analysis 3.1 One-dimensional parameter scans

Values of states/variables as a function of transition rates
We can generate lineplots of parameter scans grouped by the transition rates with the function fcn onedim parscan plot by para as shown below for the mammalian cell cycle model. This plot can be produced either for the value of model states (SI Fig. 9) or variables. On each subplot we have the stationary probability of those states or variables that show variation above a user-defined threshold as a function of the given transition rate.
With the function fcn onedim parscan plot parsensit we plot the one-dimensional parameter scan with the values of one state/variable on each subplot as a function of the transition rates, for example on Fig. 4 and 5 of the main text.

Local sensitivity of variables to transition rates
Besides the stationary value of the model's states/nodes it can be informative to plot the local sensitivity of the stationary solutions to the transition rates by visualizing their response coefficients, a metric from metabolic control analysis [6]. The response coefficient R x p of a variable x with respect to a parameter p is its derivative normalized by the local value of the variable and the parameter: As we do not compute stationary solutions symbolically we do not have the symbolic derivatives, but they can be numerically approximated by ∂x ∂p p x ≈ ∆x ∆p p x This metric tells us how a variable reacts in relative terms to a small change in a parameter. If R x p = 1 (-1) this means a 1% change in p induces a 1% increase (decrease) in x.
On SI Fig. 10 we show the response coefficients for the EMT model, highlighting stronger responses for the parameters d p63 73, u p53 and d Notch pthw consistently with Fig. 5 of the main text that shows the values of the same variables.

Multidimensional parameter sensitivity analysis with uniform resolution
With the function fcn calc paramsample table we can perform multidimensional parameter sampling with a uniform resolution. These parameter scans can be computationally expensive, eg. if the number of value for each parameter is 5, for n parameters we need 5 n calculations. For two dimensions we can visualize the results by the function fcn plot twodim parscan. Synergistic effect between transition rates can be visualized, as in SI Fig. 11 for u p53 and u Notch pthw with respect to Metastasis in the EMT model of [4].

Multidimensional parameter sensitivity analysis by Latin Hypercube Sampling (LHS)
Latin Hypercube Sampling (LHS) is a method to efficiently explore a high-dimensional parameter space. In ExaStoLog we can perform LHS by the function fcn multidim parscan latinhypcube, defining the sample size, as well as the type, mean and standard deviation of the distribution we sample from. From one-dimensional parameter scans (with the function fcn onedim parscan plot parsensit) we have the transition rates where there is significant variation in the stationary solutions and perform LHS for only these. Results of the LHS are analyzed visually and statistically by the functions described below.

Scatter plots of LHS
With the function fcn multidim parscan scatterplot we can visualize the stationary probability values (of model states or variables) from LHS as a function of transition rates as a scatterplot. The mean value of the stationary probabilities is shown by the red trend line. To calculate the mean the user needs to select the number of bins across the range of values for transition rates.

Correlations between stationary probability values of model variables
Pearson correlation coefficients between the stationary probability values of variables following LHS for the [4] model are shown on SI Fig. 13. The calculation and plotting is by the function fcn multidim parscan parvarcorrs The visually discernible trends for d BAD and u AKT are confirmed by Sobol sensitivity analysis shown on SI Fig. 15 in ExaStoLog. This metric is useful for model reduction: some variables show 100% correlation that suggests that these can be merged into single variables.

Coefficient of determination (R 2 ) between the transition rates and the stationary probability values of nodes/states
With the function fcn multidim parscan parvarcorrs we can also perform linear regression on the stationary probability values of model states/variables as a function of the transition rates. In the case of lognormal sampling for the LHS, it is recommended to do the regression as a function of the logarithm of transition rates (a built-in feature of the function). The R 2 value for each parameter-variable/state pair indicates to what extent the given transition rate predicts the stationary value. Values for the [4] model are shown on SI Fig. 14. This measure assumes a monotonic effect of transition rates on the stationary probability values, which is usually, but not necessarily the case. An alternative measure that does not have this problem is the Sobol sensitivity index. In the case of transition rates with a monotonic effect these two measures show similar, but not identical values as shown on Fig. 6 in the main text.

Sobol sensitivity index
Sobol sensitivity index is a global parameter sensitivity index that also quantifies nonlinear and nonmonotonic effects of parameters on variables. Following LHS the Sobol sensitivity index is calculated by splitting the matrix of randomly generated transition rates from the LHS into two and in one of the submatrices replacing one column (corresponding to a transition rate) with newly generated random values. Following that we recalculate the stationary solution of the model. The Sobol sensitivity index quantifies how much of the total variance in a variable is due to variation in the transition rate that was changed. The formulas for the usual numerical approximation are described in [7]. A comparison of the Sobol sensitivity index with R 2 values can be seen on Fig. 6. The result of Sobol sensitivity analysis for the [3] breast cancer

Parameter fitting
We can fit the stationary probability value of either the model's attractor states or its variables, the latter being linear combinations of the former. With the function fcn handles fitting of ExaStoLog we create anonymous functions (fcn statsol sum sq dev, fcn statsol values) to calculate the sum of squared error (SSE) of the value of the model's states/variables with respect to a vector of values (data) provided by the user. The transition rates that we want to fit are selected by the user.

Parameter fitting by simulated annealing
We use a simulated annealing algorithm anneal from MATLAB Central, with some modifications of the script for ExaStoLog (see Tutorial on GitHub). Simulated annealing is a gradient-free, probabilistic parameter fitting method. It often shows slow convergence, requiring thousands of iterations (recalculations of the stationary solution). An example is shown on Fig. 7 of the main text. Results from parameter fitting are plotted by the function fcn plot paramfitting.

Parameter fitting by numerical gradient descent
For the models we analyzed in this paper the effect of transition rates on the stationary probability value of model variables is monotonic. Therefore in some cases we can perform parameter fitting by calculating an initial gradient for the sum of squared error as a function of transition rates and reduce the SSE down to a defined level by incrementing the transition rates by this initial direction. At least in the case of some models, such as the [4] EMT model, this method is faster and leads to a larger error reduction than simulated annealing, as shown by SI Fig. compared to Fig. 7 of the main text. This method is implemented in the function fcn num grad descent, that automatically stops (after 2 steps) if the error starts to grow or stagnate. Plotting is with the same function as for simulated annealing, fcn plot paramfitting.

Relation to Probabilistic Boolean Networks (PBN)
A related but not identical approach to asynchronous-stochastic logical models (ASCTLM: asyncrhonousstochastic continuous time logical models) are probabilistic Boolean networks (PBN) [8,9]. The PBN formalism also has a stochastic aspect, as there are multiple rules per network node and one of these is stochastically selected at each updating step. Therefore a PBN is a form of ensemble modeling with one PBN model containing a number of constituent Boolean models one of which is selected to update nodes at each discrete time step. In contrast, we are calculating the stationary solution for a single model with one logical rule for each node. There are two further important differences of PBNs with respect to the ASCTLM formalism used here. First, in current implementations of PBNs [10,11,12] nodes are synchronously updated, assuming an identical timescale of activation-deactivation for the model nodes. For this reason, within the PBN formalism multiple transitions are possible from a given Boolean state, in other words there are forks in their STG with different probability weights on the outgoing arrows. However, whereas forks in the STG of PBN models arise from stochastically selecting from alternative logical rules, in the ASCTLM formalism the forks are due to individual updating of model nodes with a single logical rule for each node. We regard our ASCTLM formalism as an approximation of the different time-scales of different cellular processes, whereas in PBN the stochastic selection between multiple logical rules is more an approximation of model uncertainty when fitting to incomplete data. Second, in the current implementations of PBNs models are implemented with random perturbations [13], ie. the state (0/1) of network nodes is flipped with a given probability at any updating step. This means the system can escape attractors and makes the underlying Markov chain ergodic, ie. having a unique stationary distribution from any initial condition. We do not apply this assumption, therefore the Markov chain of our ASCTLMs is non-ergodic and the stationary solution is history-dependent. Therefore different initial conditions lead to different stationary solutions as the system cannot escape attractors. Since we are often modeling irreversible cell fate transitions, this non-ergodicity [14] is important for our purposes. To our knowledge all current implementations [9] of PBNs apply the ergodic assumption. The most scalable implementations of PBNs [11,12] use efficient stochastic simulations within the ergodic framework and can enable simulations of models with hundreds of nodes. A number of exact methods were also proposed that are based on either the power method (exponentiation of the transition matrix) [13], used for models smaller than 15 nodes, or LU decomposition of the perturbed transition matrix [10]. The latter method was scaled up to 30 nodes (2 hours of CPU time), but again was implemented using random perturbations making the entire STG strongly connected.
Due to these differences in the specification of logical rules, the updating formalism and (non-)ergodicity, the results from PBNs and ASCTLMs are not identical and the two formalisms have somewhat different scope and purpose. At the same time, the PBN formalism without ergodicity (ie. without random perturbations) is straightforward to implement in our matrix calculation framework. Namely, in this case the entry K i,j of the kinetic matrix K is the probability weight of the logical rule (applied simultaneously to all model nodes) that transforms the Boolean state i to j.