BACKGROUND
the control of gene expression is tightly regulated by complex gene regulatory networks to achieve cell type specific expression, for example in embryonic  <cit>  and blood development  <cit> .
moreover, dysregulation of gene expression can lead to disease development, including malignant disease such as leukaemia  <cit> .
a better understanding of gene regulatory networks will therefore not only advance our understanding of fundamental biological processes such as tissue development, but also provide mechanistic insights into disease processes.
the earlier versions of high-throughput expression profiling techniques were limited to measuring average gene expression across large pools of cells.
by contrast, recent technological improvements have made it possible to perform expression profiling in single cells .
protocols for the single-cell equivalent of microarray  <cit> , qpcr  <cit>  and rna sequencing  <cit>  have been developed.
one of the key advantages of single cell expression profiling is that it enables the analysis of cells that are rare in number, such as tissue stem cells.
in addition, obtaining the expression profiles of single cells is very useful for dissecting the heterogeneity within seemingly homogenous cell populations .
because single cell analysis commonly reports expression states for hundreds of individual cells, this unique information offers new opportunities for the development of algorithms that can reconstruct gene regulatory networks.
many network inference algorithms are available  <cit> , which are based on regression, correlation, mutual information and bayesian networks.
however, most of these network inference algorithms only generate a network with static representation of gene interactions.
in contrast, changes in network dynamics can be described by using dynamic models, which possess different levels of granularity and precision ranging from the simpler boolean models to more complex differential equation-based models.
more complex models such as differential equation-based models offer high precision predictions, and have been used to describe gene regulatory networks .
however, such models rely on a higher number of parameters which are often difficult to obtain and verify.
in contrast, a boolean model is one of the simplest models that can describe the dynamics of a system without the need of many parameters .
in a boolean model, each gene can take a value of  <dig> or  <dig>  which represents the absence or presence of gene expression respectively.
the interactions among genes in a boolean model are described by boolean operators like and, or and not, which closely resembles how biologists describe such interactions.
boolean models were first used to study gene regulatory networks by kauffman in the 1970s, and since then have been used extensively to study different biological systems .
while single-cell expression data offers the advantage of capturing expression profiles at single cell resolution, single-cell expression data are noisier than conventional bulk analysis.
the technical noise in single-cell expression data arises due to the low amount of input mrnas in a single cell.
this leads to two major sources of technical noise, which are pcr amplification bias and drop-outs  <cit> .
drop-outs in particular, which represent false negatives where genes are recorded as not expressed due to the low efficiency of mrna capture from single cells, represent a substantial portion of the technical noise in single-cell expression data.
therefore, network inference techniques that are robust to the effect of drop-outs are required when reconstructing networks using single-cell expression data.
boolean models are relatively robust to the presence of drop-outs due to the binarisation of expression values.
two recent studies reported algorithms for inferring boolean models from single-cell expression data  <cit> .
chen et.
al.
developed singcellnet, which uses a genetic algorithm to construct probabilistic boolean models from expected trajectories through cell states  <cit> .
however, singcellnet is restricted to small networks with less than  <dig> genes, and it only determines the network structure and transition probabilities from single-cell expression data.
the boolean rules in singcellnet are constructed via manual curation from the literature.
in another study, scns was developed by moignard et.
al.
to infer an asynchronous boolean model by analysing trajectories through a state transition graph  <cit> .
in order to infer a boolean model using scns, a connected state transition graph is required, which can be difficult to obtain from single-cell expression data.
this is because the higher the number of genes to be included in scns, the more cells will be required to build a connected state transition graph.
in addition, scns can only infer network structure by using discretised expression data, which not only leads to the loss of information, but also makes scns sensitive to the discretisation method used.
lastly, both singcellnet and scns rely on known general trajectories through the cell states, which require single-cell expression data from at least two cell types with known relationships.
here, we present a model learning algorithm btr , that is able to reconstruct and train asynchronous boolean models using single-cell expression data.
btr differs from other algorithms described above in that it can infer both network structure and boolean rules without needing information on trajectories through cell states.
we developed a scoring function based on the boolean framework, which performed favourably in comparison to a scoring function for bayesian network.
we show that btr outperforms other network inference algorithms when initial networks are supplied.
lastly, we demonstrate the capability of btr by training boolean models using single-cell qpcr and rna-seq data from haematopoietic studies.
RESULTS
a framework for scoring boolean models with single cell expression data
a boolean model b is made up of n genes x <dig>  …, xn and n update functions f <dig>  …, fn : { <dig>  1}n → { <dig>  1} each associated with a gene .
each gene can take a value x ∈ { <dig>  1}, which represents the absence or presence of gene expressions.
each update function f is expressed in terms of boolean logic by specifying the relationships among genes x <dig>  …, xn using boolean operators and , or  and not .
the main difference of asynchronous with other boolean models is the update scheme used during simulation.
an asynchronous boolean model uses the asynchronous update scheme, which specifies that at most one gene is updated between two consecutive states.
asynchronous updating is critical when modelling developmental systems that generate distinct differentiated cell types from a common progenitor, because synchronous updating generates fully deterministic models and therefore cannot capture the ability of a stem cell to mature into multiple different tissue cells.fig.
<dig> boolean model, asynchronous simulation and the framework underlying btr.
a a boolean model can be expressed graphically in terms of nodes and edges, as well as in tabular form in terms of update functions.
note that the small black node refers to and interaction.
b the asynchronous update scheme is best explained with the use of a graph representation of state space, in which each connected state differs in only one node.
starting from the initial state s
1 = { <dig>   <dig>   <dig>  1} and evaluated using the update functions in , asynchronous simulation produces a model state space with  <dig> states.
the initial state is shown in red node, while the final steady state is shown in pink node.
c the framework underlying btr.
a boolean model can be simulated to give a model state space, while a single-cell expression data can be preprocessed to give a data state space.
boolean state space scoring function can then calculate the distance score between the model and data state spaces.
lastly, btr uses the computed distance score to guide the improvement of the boolean model through an optimisation process that minimises the distance between model and data state spaces
a state in a boolean model b is represented by a boolean vector st = {x1t, …, xnt} at simulation step t. states can be generated from an initial state by systematically changing one variable at each step according to the boolean function associated with that variable.
if a state has already been encountered earlier, it is ignored.
this results in a directed graph of states as exemplified in fig.
1b, where any two connected states change in just one variable.
when all the states in the directed graph are taken together, they represent a model state space.
the initial state used in a simulation can be obtained from the expression values at time =  <dig> for a time-series expression dataset, or it can be obtained from the expression values of known parental cell types.
of note, the model state space of an asynchronous boolean model closely resembles a single-cell expression data.
the model state space contains predicted expression states that are dictated by a known gene network that underlies a boolean model; while the single-cell expression data can be viewed as a data state space which contains observed expression states that are dictated by an unknown gene network.
by fine-tuning the network rules underlying the boolean model, it should be possible to produce a predicted model state space that closely resembles an observed data state space, thereby allowing us to reconstruct the unknown gene network.
btr uses this framework to reconstruct a boolean model from single-cell expression data .
in this framework, a boolean model is represented by its model state space, while a single-cell expression dataset is represented by its data state space.
by utilising the novel boolean state space  scoring function , btr evaluates how well a particular boolean model explains the single-cell expression data by scoring the model state space with respect to the data state space.
during the model training process, btr uses a swarming hill climbing strategy to generate minimally modified boolean models based on an initial boolean model.
these minimally modified boolean models are then scored using the bss scoring function, and btr selects the best scoring boolean models for the next iteration.
by performing this process iteratively, btr reconstructs the asynchronous boolean model that can best explain a single-cell expression dataset.
boolean state space scoring represents a powerful scoring function for boolean models
how well btr performs depends heavily on the performance of the bss scoring function.
among different modelling frameworks, the bayesian network framework is known to possess several well-established scoring functions that evaluate how well a particular network fits a given dataset.
these scoring functions include log-likelihood, bayesian information criterion , bayesian dirichlet and k <dig> .
since expression data have continuous values for gene expressions, we have selected the bic scoring function, which can handle continuous variables, as a scoring function from the bayesian network framework for comparison purpose.
bss and bic scoring functions were evaluated using synthetic data.
the true network and expression data in the synthetic data were generated using genenetweaver  <cit> , which is also used in the dream <dig> network inference challenge  <cit> .
in order to simulate the zero-inflated property of single-cell expression data due to the presence of drop-outs, we introduced zero inflation into the synthetic data as described in the methods section.
an ideal scoring function should give an increasing distance score, as the evaluated network becomes increasingly different from the true network.
in order to test this, we generated a list of modified networks that are increasingly different from the true network in terms of edges.
as bayesian networks and boolean frameworks imposed different network structure constraints, the modified networks were generated separately to give a list of modified bayesian networks and another list of modified boolean networks.
although the modified bayesian and boolean networks are not identical, they possess the same number of differing edges when compared to the true network, ranging from  <dig> edges up to  <dig> differing edges.
five independent benchmark data, each with a different true network, true data and modified models, were used in the evaluation of scoring functions.
by evaluating networks using zero-inflated synthetic data, both bss and bic scoring functions performed well when acyclic networks are considered .
both scoring functions were able to give increasing distance scores as the underlying networks become increasingly different from the true network.
the bss scoring function achieves this by considering the input expression data as a data state space, and then computing the distance score by comparing the data state space with the model state space simulated from a given network.
it is expected that as a network become increasingly different, its model state space will become increasingly different from the data state space, which is reflected in the distance score as shown in fig.
2c.
to the best of our knowledge, this is the first time a scoring function that is based entirely on the boolean modelling framework has been demonstrated to give comparable performance with a scoring function for bayesian networks.fig.
<dig> bss scoring function compares favourably with bic scoring function on acyclic networks.
a acyclic networks generated from genenetweaver that are designated as the true acyclic networks.
each node corresponds to a gene.
black edges indicate activation interactions, while red edges indicate inhibition interactions.
mean distance scores computed using b bic scoring function and c bss scoring function for modified networks that are increasingly different from the true network in terms of edges using zero-inflated synthetic expression data.
the modified networks contain from two edges up to forty different edges when compared with the true network.
each data point is the mean distance score of  <dig> different random modified networks that contain the same number of different edges with respect to the true network.
the error bar is the standard error of the mean
as indicated in the results for network  <dig> , the bss scoring function is dependent on the underlying true network structure in certain cases and will work better on distinguishing networks that are very different.
however the bss scoring function has a distinct advantage over scoring functions for bayesian networks.
the bayesian networks are known to impose relatively strict constraints on permissible network structures, in particular bayesian networks are not allowed to contain any cyclic network structure.
therefore scoring functions for bayesian networks cannot be used to evaluate cyclic networks.
cyclic networks are ubiquitous in biological systems, in which cyclic motifs can be present in the form of negative and positive feedback loops.
boolean models on the other hand are allowed to have any number of cyclic motifs in the networks.
therefore, the bss scoring function can be used to compute scores for cyclic networks.
by using another five independent benchmark data with true networks that contain at least one cycle, the distance scores for modified networks were computed .
the distance scores for cyclic networks have more fluctuations compared to acyclic networks due to the presence of cyclic motifs.
however, the general trend where the distance scores increase as the underlying networks become increasingly different from the true network was still observed.fig.
<dig> bss scoring function is able to calculate distance scores for cyclic networks.
a cyclic networks generated from genenetweaver that are designated as the true cyclic networks.
each node corresponds to a gene.
black edges indicate activation interactions, while red edges indicate inhibition interactions.
b mean distance scores computed using bss scoring function for modified networks that are increasingly different from the true network in terms of edges using zero-inflated synthetic expression data.
the modified networks contain from two edges up to forty different edges when compared with the true network.
each data point is the mean distance score of  <dig> different random modified networks that contain the same number of different edges with respect to the true network.
the error bar is the standard error of the mean
we have also evaluated the series of acyclic and cyclic networks using non zero-inflated data .
when the results computed with non zero-inflated data are compared to the results computed using zero-inflated data, we can see that zero-inflation has no effect on bic scores and a small effect on bss scores that does not affect the general trend .
in summary, the relative mean scores that average across the results of all networks  show that although the bic scoring function performs slightly better than the bss scoring function, the bss scoring function has the advantage that it can evaluate cyclic networks.fig.
<dig> summary of bic and bss scoring functions.
mean scores have been calculated across all networks  for bic and bss scoring functions calculated using zero-inflated synthetic expression data.
all scores have been standardised for comparison purpose, such that the scores range from  <dig> to 1
btr accurately infers the networks underlying synthetic datasets
next, we compared the network inference performance of btr with other well-known network inference algorithms.
two search algorithms guided by the bss boolean and bic bayesian network scoring functions were included in the comparison, indicated as btr and bic respectively.
the search algorithms used for both scoring functions are both based on hill climbing.
the additional network inference algorithms included in the comparison are bestfit  <cit> , aracne  <cit> , clr  <cit> , bc3net  <cit> , genenet  <cit>  and genie <dig>  <cit>  .
by using the same synthetic networks, as well as both non zero-inflated and zero-inflated synthetic data, we performed network inference using the synthetic expression data alone without any extra information.
in contrast to the dream <dig> challenge  <cit>  which also provides perturbed expression data, only a single type of expression data is provided to all the network inference algorithms, which is the wild type time course expression data in steady state.
for btr, besides performing inference with only expression data , we also performed inference with both expression data and initial networks  to show that btr is able to use initial networks with known network structure to improve the inference process.
the initial networks are generated randomly to contain  <dig> edges that are different compared with the true networks.
the performance of the network inference algorithms is assessed in terms of f-scores  <cit>  .
in order to allow comparisons on the performance across all network inference algorithms tested, we calculated the f-scores based only on the presence or absence of edges, while ignoring any additional information such as the types of edges.fig.
<dig> btr outperforms other network inference algorithms.
mean f-scores of network inference algorithms inferred using a non zero-inflated synthetic data and b zero-inflated synthetic data.
ten true synthetic networks  were used in the assessment of these network inference algorithms.
plots titled ‘both’ show the combined results of acyclic and cyclic network inference.
the error bar is the standard error of the mean
in terms of acyclic networks, the results show that the top inference algorithms using either non zero-inflated or zero-inflated data are btr-wi, clr, bic and btr-wo.
as for cyclic networks, the top inference algorithms differ between using non zero-inflated and zero-inflated data.
btr-wi, btr-wo, clr and bc3net gave the best performance with non zero-inflated data, while btr-wi, aracne, genie <dig> and clr gave the best performance with zero-inflated data.
when all results are taken together, btr-wi, clr, btr-wo and genie <dig> gave the best performance overall.
note that the ranking of network inference algorithms in this study differs from the ranking of the dream study because different scoring criteria are used  and receiver operating characteristic  curves in the dream study); and the dream study was done using multiple types of synthetic data, such as expression data with gene perturbations.
in general, the presence of drop-outs affects the performance of network inference algorithms in different ways .
in cases such as bc3net and genenet, their performance decreases when drop-outs are present, while the impact of drop-outs on the performance of btr is minimal.
interestingly, the performance of bestfit increases with the presence of drop-outs, possibly due to better binarisation of data due to the information given by drop-outs.
as both btr and bestfit are algorithms for inferring boolean model, this result provides further support that boolean models are robust to the presence of drop-outs in single-cell expression data.
when given an initial network as in btr-wi, the btr algorithm was able to perform very well in locating the true network.
while the performance of the btr algorithm without an initial network  is comparable with other inference algorithms, btr-wo scored less well compared to btr-wi.
this indicates that the greedy hill climbing search strategy implemented in btr may not be able to traverse the solution space efficiently without any initial information.
taken together, while btr can be used for reconstructing network models without initial information, btr performed the best when it is used to train and improve on existing networks that contain a partially true structure.
it is also worth noting that btr produced a dynamic model with a directed underlying static network, in contrast to most other algorithms such as clr that only produce an undirected static network.
btr predicts gene interactions by training haematopoietic boolean models
we next wanted to apply btr to biological data to evaluate its utility to biologists.
haematopoiesis research has provided many paradigms for modern biological research, and was one of the first fields to embrace single cell expression profiling  <cit> .
moreover, literature curated boolean network models have been reported both for blood stem cell maintenance and blood progenitor differentiation  <cit> .
the single-cell expression data used here includes single-cell qpcr and single-cell rna-seq data, which are both obtained from  <cit> .
the two boolean models will be referred to as the bonzanni model  <cit>   and the krumsiek model  <cit>  .
both models had been constructed via manual literature curation by the authors of the original papers.
the bonzanni model aimed to capture haematopoietic stem cell  self-renewal capacity, while the krumsiek model describes the differentiation process of the erythro-myeloid lineage in haematopoiesis.fig.
<dig> btr predicts gene interactions by training the bonzanni and krumsiek boolean models.
a original bonzanni model.
b trained bonzanni model.
c original krumsiek model.
d trained krumsiek model.
round orange nodes indicate genes, square black nodes indicate and gates that combine the two input gene interactions.
blue edges indicate activation interactions, red edges indicate inhibition interactions.
dashed lines in the original models indicate edges that are present in the original models, but are removed in the trained models.
dashed lines in the trained models indicate edges that are added to the trained models and are not present in the original models
we firstly trained the bonzanni model using single-cell rna-seq data collected from hscs.
compared to the original model, the resulting trained bonzanni model  shows the deletions of ten gene interactions and the additions of thirteen gene interactions .
the state space of the trained bonzanni model contains  <dig> states when simulated using the initial state used in the original study .
of note, there are many densely connected transitional states in the state space, which may be related to the complexity of cell fate decision making processes in multipotent progenitor cells.
steady state analysis performed showed that the steady states of the trained bonzanni model are almost identical to the steady states of the original bonzanni model , except with the absence of cyclic steady states.
the authors suggested that the cyclic steady states in the original bonzanni model correspond to the self-renewal maintenance loop in hscs, which is not present in our trained model possibly because the number of cells profiled by single-cell rna-seq is not enough to sufficiently capture the hsc self-renewal expression signature.
we then trained the krumsiek model by using single-cell qpcr data collected from over  <dig> cells along the erythro-myeloid lineage, which includes common myeloid progenitors, granulocyte-monocyte progenitors and myeloid-erythroid progenitors.
in order to demonstrate that btr can be used in cases where we may want to extend a current boolean model by adding more genes to it, we have used btr to train and add two additional genes to the krumsiek model.
the resulting trained krumsiek model  contains three deleted gene interaction and twelve added gene interactions  when compared to the original krumsiek model.
for the two additional genes ldb <dig> and lmo <dig>  btr has predicted gene interactions among ldb <dig>  lmo <dig>  fli <dig>  gata <dig> and gata <dig>  previous studies have shown that genome-wide binding profiles for lmo <dig>  gata <dig> and fli <dig> show significant overlaps  <cit> , and that ldb <dig> also occupies nearly all of the binding sites of gata <dig>  <cit> , consistent with a model where these tfs engage in combinatorial interactions.
the state space of the trained krumsiek model contains  <dig> states when simulated using the initial state used in the original study .
the two steady states reachable in this state space may correspond well to cell populations that are primed for the erythrocyte and myeloid lineage divergence.
when examining the steady states reachable from all possible initial states, the trained krumsiek model produces additional steady states when compared with the original model due to the addition of two extra genes , which may correspond to intermediate cell types along the erythro-myeloid differentiation pathway.table  <dig> differences in gene interactions between original and trained bonzanni and krumsiek models
gata <dig> and hhex inhibits gata2
scl and gata <dig> activates scl
fli <dig> and gata <dig> activates scl
gata <dig> and scl activates fli1
fli <dig> activates runx1
fli <dig> activates erg
erg activates erg
gata <dig> and fog <dig> inhibits gata2
sfpi <dig> and gata <dig> inhibits gata1
sfpi <dig> activates sfpi1
sfpi <dig> and erg activates sfpi1
hhex and runx <dig> inhibits gata2
eto <dig> and hhex inhibits gata2
sfpi <dig> activates hhex
gata <dig> inhibits gata2
fog <dig> inhibits smad6
fog <dig> inhibits fli1
eto <dig> activates gata1
eto <dig> activates fog1
hhex inhibits gata1
hhex activates eto2
hhex inhibits erg
fli <dig> inhibits runx1
cebpa activates gfi1
gata <dig> inhibits sfpi1
sfpi <dig> inhibits gata2
gfi <dig> activates fli1
cebpa activates fli1
fli <dig> and gata <dig> inhibits sfpi1
fli <dig> and sfpi <dig> inhibits gata2
cebpa inhibits gata1
egrnab activates ldb1
cjun activates ldb1
ldb <dig> inhibits gata1
gata <dig> activates ldb1
ldb <dig> activates lmo2
ldb <dig> and gata <dig> activates lmo2
lmo <dig> inhibits gata1
other gene interactions that were not modified are not listed in this table.
each gene interaction corresponds to an edge on the network
fig.
<dig> state spaces for the trained bonzanni and krumsiek boolean models.
a state space of trained bonzanni model.
b state space of trained krumsiek model.
blue nodes represent transitional model states, while pink nodes represent steady model states.
each arrow indicates transitions among states
fig.
<dig> steady states for the bonzanni and krumsiek boolean models.
a steady states of bonzanni models.
both original and trained models contain two point steady states and one cyclic steady state each.
b steady states of krumsiek models.
original model contains six point steady states, while trained model contains nine point steady states.
black box indicates expression is absent , red box indicates expression is present , purple box indicates expression can be absent or present .
white box is used to indicate that the additional genes included in trained krumsiek model that are not present in the original krumsiek model
taken together, the result suggests that both the trained bonzanni and krumsiek models have been trained by btr to predict new gene interactions which give rise to interesting state spaces and steady state properties.
note that the state space of the trained bonzanni model is substantially larger than the state space of the trained krumsiek model due to the denser interactions among genes and a lower proportion of inhibitory edges in the trained bonzanni model .
CONCLUSIONS
we have developed the btr model learning algorithm for training asynchronous boolean models using single-cell expression data.
the key component in btr is a novel boolean state space  scoring function, which btr uses to infer a boolean model through an optimisation process.
we have shown that the new bss scoring function is capable of giving meaningful scores to networks when compared with the bic scoring function for bayesian networks.
we then showed that when compared to other network reconstruction algorithms, btr gave the best result when initial networks were provided.
in two case studies, we have demonstrated that btr is capable of suggesting modifications to existing boolean models based on information from single-cell qpcr and rna-seq data.
finally, we anticipate btr to be a useful addition to the current toolbox for processing and understanding single-cell expression data, as it provides significant new capabilities for regulatory network modelling in a user-friendly way.

 the control of gene expression is tightly regulated by complex gene regulatory networks to achieve cell type specific expression, for example in embryonic and blood development. moreover, dysregulation of gene expression can lead to disease development, including malignant disease such as leukaemia. a better understanding of gene regulatory networks will therefore not only advance our understanding of fundamental biological processes such as tissue development, but also provide mechanistic insights into disease processes. the earlier versions of high-throughput expression profiling techniques were limited to measuring average gene expression across large pools of cells. by contrast, recent technological improvements have made it possible to perform expression profiling in single cells. protocols for the single-cell equivalent of microarray, qpcr and rna sequencing have been developed. one of the key advantages of single cell expression profiling is that it enables the analysis of cells that are rare in number, such as tissue stem cells. in addition, obtaining the expression profiles of single cells is very useful for dissecting the heterogeneity within seemingly homogenous cell populations. because single cell analysis commonly reports expression states for hundreds of individual cells, this unique information offers new opportunities for the development of algorithms that can reconstruct gene regulatory networks. many network inference algorithms are available, which are based on regression, correlation, mutual information and bayesian networks. however, most of these network inference algorithms only generate a network with static representation of gene interactions. in contrast, changes in network dynamics can be described by using dynamic models, which possess different levels of granularity and precision ranging from the simpler boolean models to more complex differential equation-based models. more complex models such as differential equation-based models offer high precision predictions, and have been used to describe gene regulatory networks. however, such models rely on a higher number of parameters which are often difficult to obtain and verify. in contrast, a boolean model is one of the simplest models that can describe the dynamics of a system without the need of many parameters. in a boolean model, each gene can take a value of or which represents the absence or presence of gene expression respectively. the interactions among genes in a boolean model are described by boolean operators like and, or and not, which closely resembles how biologists describe such interactions. boolean models were first used to study gene regulatory networks by kauffman in the 1970s, and since then have been used extensively to study different biological systems. while single-cell expression data offers the advantage of capturing expression profiles at single cell resolution, single-cell expression data are noisier than conventional bulk analysis. the technical noise in single-cell expression data arises due to the low amount of input mrnas in a single cell. this leads to two major sources of technical noise, which are pcr amplification bias and drop-outs. drop-outs in particular, which represent false negatives where genes are recorded as not expressed due to the low efficiency of mrna capture from single cells, represent a substantial portion of the technical noise in single-cell expression data. therefore, network inference techniques that are robust to the effect of drop-outs are required when reconstructing networks using single-cell expression data. boolean models are relatively robust to the presence of drop-outs due to the binarisation of expression values. two recent studies reported algorithms for inferring boolean models from single-cell expression data. chen et. al. developed singcellnet, which uses a genetic algorithm to construct probabilistic boolean models from expected trajectories through cell states. however, singcellnet is restricted to small networks with less than genes, and it only determines the network structure and transition probabilities from single-cell expression data. the boolean rules in singcellnet are constructed via manual curation from the literature. in another study, scns was developed by moignard et. al. to infer an asynchronous boolean model by analysing trajectories through a state transition graph. in order to infer a boolean model using scns, a connected state transition graph is required, which can be difficult to obtain from single-cell expression data. this is because the higher the number of genes to be included in scns, the more cells will be required to build a connected state transition graph. in addition, scns can only infer network structure by using discretised expression data, which not only leads to the loss of information, but also makes scns sensitive to the discretisation method used. lastly, both singcellnet and scns rely on known general trajectories through the cell states, which require single-cell expression data from at least two cell types with known relationships. here, we present a model learning algorithm btr, that is able to reconstruct and train asynchronous boolean models using single-cell expression data. btr differs from other algorithms described above in that it can infer both network structure and boolean rules without needing information on trajectories through cell states. we developed a scoring function based on the boolean framework, which performed favourably in comparison to a scoring function for bayesian network. we show that btr outperforms other network inference algorithms when initial networks are supplied. lastly, we demonstrate the capability of btr by training boolean models using single-cell qpcr and rna-seq data from haematopoietic studies.  a framework for scoring boolean models with single cell expression data a boolean model b is made up of n genes x  …, xn and n update functions f  …, fn: {  1}n → {  1} each associated with a gene. each gene can take a value x ∈ {  1}, which represents the absence or presence of gene expressions. each update function f is expressed in terms of boolean logic by specifying the relationships among genes x  …, xn using boolean operators and, or and not. the main difference of asynchronous with other boolean models is the update scheme used during simulation. an asynchronous boolean model uses the asynchronous update scheme, which specifies that at most one gene is updated between two consecutive states. asynchronous updating is critical when modelling developmental systems that generate distinct differentiated cell types from a common progenitor, because synchronous updating generates fully deterministic models and therefore cannot capture the ability of a stem cell to mature into multiple different tissue cells.fig. boolean model, asynchronous simulation and the framework underlying btr. a a boolean model can be expressed graphically in terms of nodes and edges, as well as in tabular form in terms of update functions. note that the small black node refers to and interaction. b the asynchronous update scheme is best explained with the use of a graph representation of state space, in which each connected state differs in only one node. starting from the initial state s 1 = {      1} and evaluated using the update functions in, asynchronous simulation produces a model state space with states. the initial state is shown in red node, while the final steady state is shown in pink node. c the framework underlying btr. a boolean model can be simulated to give a model state space, while a single-cell expression data can be preprocessed to give a data state space. boolean state space scoring function can then calculate the distance score between the model and data state spaces. lastly, btr uses the computed distance score to guide the improvement of the boolean model through an optimisation process that minimises the distance between model and data state spaces a state in a boolean model b is represented by a boolean vector st = {x1t, …, xnt} at simulation step t. states can be generated from an initial state by systematically changing one variable at each step according to the boolean function associated with that variable. if a state has already been encountered earlier, it is ignored. this results in a directed graph of states as exemplified in fig. 1b, where any two connected states change in just one variable. when all the states in the directed graph are taken together, they represent a model state space. the initial state used in a simulation can be obtained from the expression values at time =  for a time-series expression dataset, or it can be obtained from the expression values of known parental cell types. of note, the model state space of an asynchronous boolean model closely resembles a single-cell expression data. the model state space contains predicted expression states that are dictated by a known gene network that underlies a boolean model; while the single-cell expression data can be viewed as a data state space which contains observed expression states that are dictated by an unknown gene network. by fine-tuning the network rules underlying the boolean model, it should be possible to produce a predicted model state space that closely resembles an observed data state space, thereby allowing us to reconstruct the unknown gene network. btr uses this framework to reconstruct a boolean model from single-cell expression data. in this framework, a boolean model is represented by its model state space, while a single-cell expression dataset is represented by its data state space. by utilising the novel boolean state space scoring function, btr evaluates how well a particular boolean model explains the single-cell expression data by scoring the model state space with respect to the data state space. during the model training process, btr uses a swarming hill climbing strategy to generate minimally modified boolean models based on an initial boolean model. these minimally modified boolean models are then scored using the bss scoring function, and btr selects the best scoring boolean models for the next iteration. by performing this process iteratively, btr reconstructs the asynchronous boolean model that can best explain a single-cell expression dataset. boolean state space scoring represents a powerful scoring function for boolean models how well btr performs depends heavily on the performance of the bss scoring function. among different modelling frameworks, the bayesian network framework is known to possess several well-established scoring functions that evaluate how well a particular network fits a given dataset. these scoring functions include log-likelihood, bayesian information criterion, bayesian dirichlet and k. since expression data have continuous values for gene expressions, we have selected the bic scoring function, which can handle continuous variables, as a scoring function from the bayesian network framework for comparison purpose. bss and bic scoring functions were evaluated using synthetic data. the true network and expression data in the synthetic data were generated using genenetweaver, which is also used in the dream network inference challenge. in order to simulate the zero-inflated property of single-cell expression data due to the presence of drop-outs, we introduced zero inflation into the synthetic data as described in the methods section. an ideal scoring function should give an increasing distance score, as the evaluated network becomes increasingly different from the true network. in order to test this, we generated a list of modified networks that are increasingly different from the true network in terms of edges. as bayesian networks and boolean frameworks imposed different network structure constraints, the modified networks were generated separately to give a list of modified bayesian networks and another list of modified boolean networks. although the modified bayesian and boolean networks are not identical, they possess the same number of differing edges when compared to the true network, ranging from edges up to differing edges. five independent benchmark data, each with a different true network, true data and modified models, were used in the evaluation of scoring functions. by evaluating networks using zero-inflated synthetic data, both bss and bic scoring functions performed well when acyclic networks are considered. both scoring functions were able to give increasing distance scores as the underlying networks become increasingly different from the true network. the bss scoring function achieves this by considering the input expression data as a data state space, and then computing the distance score by comparing the data state space with the model state space simulated from a given network. it is expected that as a network become increasingly different, its model state space will become increasingly different from the data state space, which is reflected in the distance score as shown in fig. 2c. to the best of our knowledge, this is the first time a scoring function that is based entirely on the boolean modelling framework has been demonstrated to give comparable performance with a scoring function for bayesian networks.fig. bss scoring function compares favourably with bic scoring function on acyclic networks. a acyclic networks generated from genenetweaver that are designated as the true acyclic networks. each node corresponds to a gene. black edges indicate activation interactions, while red edges indicate inhibition interactions. mean distance scores computed using b bic scoring function and c bss scoring function for modified networks that are increasingly different from the true network in terms of edges using zero-inflated synthetic expression data. the modified networks contain from two edges up to forty different edges when compared with the true network. each data point is the mean distance score of different random modified networks that contain the same number of different edges with respect to the true network. the error bar is the standard error of the mean as indicated in the results for network, the bss scoring function is dependent on the underlying true network structure in certain cases and will work better on distinguishing networks that are very different. however the bss scoring function has a distinct advantage over scoring functions for bayesian networks. the bayesian networks are known to impose relatively strict constraints on permissible network structures, in particular bayesian networks are not allowed to contain any cyclic network structure. therefore scoring functions for bayesian networks cannot be used to evaluate cyclic networks. cyclic networks are ubiquitous in biological systems, in which cyclic motifs can be present in the form of negative and positive feedback loops. boolean models on the other hand are allowed to have any number of cyclic motifs in the networks. therefore, the bss scoring function can be used to compute scores for cyclic networks. by using another five independent benchmark data with true networks that contain at least one cycle, the distance scores for modified networks were computed. the distance scores for cyclic networks have more fluctuations compared to acyclic networks due to the presence of cyclic motifs. however, the general trend where the distance scores increase as the underlying networks become increasingly different from the true network was still observed.fig. bss scoring function is able to calculate distance scores for cyclic networks. a cyclic networks generated from genenetweaver that are designated as the true cyclic networks. each node corresponds to a gene. black edges indicate activation interactions, while red edges indicate inhibition interactions. b mean distance scores computed using bss scoring function for modified networks that are increasingly different from the true network in terms of edges using zero-inflated synthetic expression data. the modified networks contain from two edges up to forty different edges when compared with the true network. each data point is the mean distance score of different random modified networks that contain the same number of different edges with respect to the true network. the error bar is the standard error of the mean we have also evaluated the series of acyclic and cyclic networks using non zero-inflated data. when the results computed with non zero-inflated data are compared to the results computed using zero-inflated data, we can see that zero-inflation has no effect on bic scores and a small effect on bss scores that does not affect the general trend. in summary, the relative mean scores that average across the results of all networks show that although the bic scoring function performs slightly better than the bss scoring function, the bss scoring function has the advantage that it can evaluate cyclic networks.fig. summary of bic and bss scoring functions. mean scores have been calculated across all networks for bic and bss scoring functions calculated using zero-inflated synthetic expression data. all scores have been standardised for comparison purpose, such that the scores range from to 1 btr accurately infers the networks underlying synthetic datasets next, we compared the network inference performance of btr with other well-known network inference algorithms. two search algorithms guided by the bss boolean and bic bayesian network scoring functions were included in the comparison, indicated as btr and bic respectively. the search algorithms used for both scoring functions are both based on hill climbing. the additional network inference algorithms included in the comparison are bestfit, aracne, clr, bc3net, genenet and genie. by using the same synthetic networks, as well as both non zero-inflated and zero-inflated synthetic data, we performed network inference using the synthetic expression data alone without any extra information. in contrast to the dream challenge which also provides perturbed expression data, only a single type of expression data is provided to all the network inference algorithms, which is the wild type time course expression data in steady state. for btr, besides performing inference with only expression data, we also performed inference with both expression data and initial networks to show that btr is able to use initial networks with known network structure to improve the inference process. the initial networks are generated randomly to contain edges that are different compared with the true networks. the performance of the network inference algorithms is assessed in terms of f-scores. in order to allow comparisons on the performance across all network inference algorithms tested, we calculated the f-scores based only on the presence or absence of edges, while ignoring any additional information such as the types of edges.fig. btr outperforms other network inference algorithms. mean f-scores of network inference algorithms inferred using a non zero-inflated synthetic data and b zero-inflated synthetic data. ten true synthetic networks were used in the assessment of these network inference algorithms. plots titled ‘both’ show the combined results of acyclic and cyclic network inference. the error bar is the standard error of the mean in terms of acyclic networks, the results show that the top inference algorithms using either non zero-inflated or zero-inflated data are btr-wi, clr, bic and btr-wo. as for cyclic networks, the top inference algorithms differ between using non zero-inflated and zero-inflated data. btr-wi, btr-wo, clr and bc3net gave the best performance with non zero-inflated data, while btr-wi, aracne, genie and clr gave the best performance with zero-inflated data. when all results are taken together, btr-wi, clr, btr-wo and genie gave the best performance overall. note that the ranking of network inference algorithms in this study differs from the ranking of the dream study because different scoring criteria are used and receiver operating characteristic curves in the dream study); and the dream study was done using multiple types of synthetic data, such as expression data with gene perturbations. in general, the presence of drop-outs affects the performance of network inference algorithms in different ways. in cases such as bc3net and genenet, their performance decreases when drop-outs are present, while the impact of drop-outs on the performance of btr is minimal. interestingly, the performance of bestfit increases with the presence of drop-outs, possibly due to better binarisation of data due to the information given by drop-outs. as both btr and bestfit are algorithms for inferring boolean model, this result provides further support that boolean models are robust to the presence of drop-outs in single-cell expression data. when given an initial network as in btr-wi, the btr algorithm was able to perform very well in locating the true network. while the performance of the btr algorithm without an initial network is comparable with other inference algorithms, btr-wo scored less well compared to btr-wi. this indicates that the greedy hill climbing search strategy implemented in btr may not be able to traverse the solution space efficiently without any initial information. taken together, while btr can be used for reconstructing network models without initial information, btr performed the best when it is used to train and improve on existing networks that contain a partially true structure. it is also worth noting that btr produced a dynamic model with a directed underlying static network, in contrast to most other algorithms such as clr that only produce an undirected static network. btr predicts gene interactions by training haematopoietic boolean models we next wanted to apply btr to biological data to evaluate its utility to biologists. haematopoiesis research has provided many paradigms for modern biological research, and was one of the first fields to embrace single cell expression profiling. moreover, literature curated boolean network models have been reported both for blood stem cell maintenance and blood progenitor differentiation. the single-cell expression data used here includes single-cell qpcr and single-cell rna-seq data, which are both obtained from. the two boolean models will be referred to as the bonzanni model and the krumsiek model. both models had been constructed via manual literature curation by the authors of the original papers. the bonzanni model aimed to capture haematopoietic stem cell self-renewal capacity, while the krumsiek model describes the differentiation process of the erythro-myeloid lineage in haematopoiesis.fig. btr predicts gene interactions by training the bonzanni and krumsiek boolean models. a original bonzanni model. b trained bonzanni model. c original krumsiek model. d trained krumsiek model. round orange nodes indicate genes, square black nodes indicate and gates that combine the two input gene interactions. blue edges indicate activation interactions, red edges indicate inhibition interactions. dashed lines in the original models indicate edges that are present in the original models, but are removed in the trained models. dashed lines in the trained models indicate edges that are added to the trained models and are not present in the original models we firstly trained the bonzanni model using single-cell rna-seq data collected from hscs. compared to the original model, the resulting trained bonzanni model shows the deletions of ten gene interactions and the additions of thirteen gene interactions. the state space of the trained bonzanni model contains states when simulated using the initial state used in the original study. of note, there are many densely connected transitional states in the state space, which may be related to the complexity of cell fate decision making processes in multipotent progenitor cells. steady state analysis performed showed that the steady states of the trained bonzanni model are almost identical to the steady states of the original bonzanni model, except with the absence of cyclic steady states. the authors suggested that the cyclic steady states in the original bonzanni model correspond to the self-renewal maintenance loop in hscs, which is not present in our trained model possibly because the number of cells profiled by single-cell rna-seq is not enough to sufficiently capture the hsc self-renewal expression signature. we then trained the krumsiek model by using single-cell qpcr data collected from over cells along the erythro-myeloid lineage, which includes common myeloid progenitors, granulocyte-monocyte progenitors and myeloid-erythroid progenitors. in order to demonstrate that btr can be used in cases where we may want to extend a current boolean model by adding more genes to it, we have used btr to train and add two additional genes to the krumsiek model. the resulting trained krumsiek model contains three deleted gene interaction and twelve added gene interactions when compared to the original krumsiek model. for the two additional genes ldb and lmo btr has predicted gene interactions among ldb lmo fli gata and gata previous studies have shown that genome-wide binding profiles for lmo gata and fli show significant overlaps, and that ldb also occupies nearly all of the binding sites of gata, consistent with a model where these tfs engage in combinatorial interactions. the state space of the trained krumsiek model contains states when simulated using the initial state used in the original study. the two steady states reachable in this state space may correspond well to cell populations that are primed for the erythrocyte and myeloid lineage divergence. when examining the steady states reachable from all possible initial states, the trained krumsiek model produces additional steady states when compared with the original model due to the addition of two extra genes, which may correspond to intermediate cell types along the erythro-myeloid differentiation pathway.table differences in gene interactions between original and trained bonzanni and krumsiek models gata and hhex inhibits gata2 scl and gata activates scl fli and gata activates scl gata and scl activates fli1 fli activates runx1 fli activates erg erg activates erg gata and fog inhibits gata2 sfpi and gata inhibits gata1 sfpi activates sfpi1 sfpi and erg activates sfpi1 hhex and runx inhibits gata2 eto and hhex inhibits gata2 sfpi activates hhex gata inhibits gata2 fog inhibits smad6 fog inhibits fli1 eto activates gata1 eto activates fog1 hhex inhibits gata1 hhex activates eto2 hhex inhibits erg fli inhibits runx1 cebpa activates gfi1 gata inhibits sfpi1 sfpi inhibits gata2 gfi activates fli1 cebpa activates fli1 fli and gata inhibits sfpi1 fli and sfpi inhibits gata2 cebpa inhibits gata1 egrnab activates ldb1 cjun activates ldb1 ldb inhibits gata1 gata activates ldb1 ldb activates lmo2 ldb and gata activates lmo2 lmo inhibits gata1 other gene interactions that were not modified are not listed in this table. each gene interaction corresponds to an edge on the network fig. state spaces for the trained bonzanni and krumsiek boolean models. a state space of trained bonzanni model. b state space of trained krumsiek model. blue nodes represent transitional model states, while pink nodes represent steady model states. each arrow indicates transitions among states fig. steady states for the bonzanni and krumsiek boolean models. a steady states of bonzanni models. both original and trained models contain two point steady states and one cyclic steady state each. b steady states of krumsiek models. original model contains six point steady states, while trained model contains nine point steady states. black box indicates expression is absent, red box indicates expression is present, purple box indicates expression can be absent or present. white box is used to indicate that the additional genes included in trained krumsiek model that are not present in the original krumsiek model taken together, the result suggests that both the trained bonzanni and krumsiek models have been trained by btr to predict new gene interactions which give rise to interesting state spaces and steady state properties. note that the state space of the trained bonzanni model is substantially larger than the state space of the trained krumsiek model due to the denser interactions among genes and a lower proportion of inhibitory edges in the trained bonzanni model.  we have developed the btr model learning algorithm for training asynchronous boolean models using single-cell expression data. the key component in btr is a novel boolean state space scoring function, which btr uses to infer a boolean model through an optimisation process. we have shown that the new bss scoring function is capable of giving meaningful scores to networks when compared with the bic scoring function for bayesian networks. we then showed that when compared to other network reconstruction algorithms, btr gave the best result when initial networks were provided. in two case studies, we have demonstrated that btr is capable of suggesting modifications to existing boolean models based on information from single-cell qpcr and rna-seq data. finally, we anticipate btr to be a useful addition to the current toolbox for processing and understanding single-cell expression data, as it provides significant new capabilities for regulatory network modelling in a user-friendly way. 
here we introduce btr, an algorithm for training asynchronous boolean models with single-cell expression data using a novel boolean state space scoring function. btr is capable of refining existing boolean models and reconstructing new boolean models by improving the match between model prediction and expression data. lastly, we introduce two case studies, in which we use btr to improve published boolean models in order to generate potentially new biological insights. btr provides a novel way to refine or reconstruct boolean models using single-cell expression data. boolean model is particularly useful for network reconstruction using single-cell data because it is more robust to the effect of drop-outs. given the simplicity of boolean models and the rapid adoption of single-cell genomics by biologists, btr has the potential to make an impact across many fields of biomedical research. asynchronous boolean modelsingle-cell gene expressionmodel learningnetwork reconstructionboolean scoring functionexecutable modelbloodwise12029http://dx.doi.org/ /501100000268biotechnology and biological sciences research councilbb/i00050x/1http://dx.doi.org/ /501100000580cambridge institute for medical research, university of cambridgerp-pg-0310-1002http://dx.doi.org/ /100006112microsoft research2012-023http://dx.doi.org/ /100004440wellcome trust100140/z/12/z; 097922/z/11/zissue-copyright-statement© the author 2016 while many algorithms exist to infer a gene regulatory network, very few of them are able to harness the extra expression states present in single-cell expression data without getting adversely affected by the substantial technical noise present. the use of single-cell expression profiling technique allows the profiling of the expression states of hundreds of cells, but these expression states are typically noisier due to the presence of technical artefacts such as drop-outs. in addition, btr does not assume any relationship in the expression states among cells, it is useful for reconstructing a gene regulatory network with as few assumptions as possible. in addition, we show that btr outperforms many other network inference algorithms in both bulk and single-cell synthetic expression data. we demonstrate that the boolean scoring function performed favourably against the bic scoring function for bayesian networks.
 rapid technological innovation for the generation of single-cell genomics data presents new challenges and opportunities for bioinformatics analysis. one such area lies in the development of new ways to train gene regulatory networks. the use of single-cell expression profiling technique allows the profiling of the expression states of hundreds of cells, but these expression states are typically noisier due to the presence of technical artefacts such as drop-outs. while many algorithms exist to infer a gene regulatory network, very few of them are able to harness the extra expression states present in single-cell expression data without getting adversely affected by the substantial technical noise present.  here we introduce btr, an algorithm for training asynchronous boolean models with single-cell expression data using a novel boolean state space scoring function. btr is capable of refining existing boolean models and reconstructing new boolean models by improving the match between model prediction and expression data. we demonstrate that the boolean scoring function performed favourably against the bic scoring function for bayesian networks. in addition, we show that btr outperforms many other network inference algorithms in both bulk and single-cell synthetic expression data. lastly, we introduce two case studies, in which we use btr to improve published boolean models in order to generate potentially new biological insights.  btr provides a novel way to refine or reconstruct boolean models using single-cell expression data. boolean model is particularly useful for network reconstruction using single-cell data because it is more robust to the effect of drop-outs. in addition, btr does not assume any relationship in the expression states among cells, it is useful for reconstructing a gene regulatory network with as few assumptions as possible. given the simplicity of boolean models and the rapid adoption of single-cell genomics by biologists, btr has the potential to make an impact across many fields of biomedical research. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords asynchronous boolean modelsingle-cell gene expressionmodel learningnetwork reconstructionboolean scoring functionexecutable modelbloodwise12029http://dx.doi.org/ /501100000268biotechnology and biological sciences research councilbb/i00050x/1http://dx.doi.org/ /501100000580cambridge institute for medical research, university of cambridgerp-pg-0310-1002http://dx.doi.org/ /100006112microsoft research2012-023http://dx.doi.org/ /100004440wellcome trust100140/z/12/z; 097922/z/11/zissue-copyright-statement© the author 2016 