Computer models are widely used to predict behaviors of biological systems, generally by simulating a number of instances of a model and enumerating the observed outcomes. Comprehensive simulations of single cells may be realistic given recent progress in constructing elaborate cellular models [1, 2], but simulation of an entire tissue is far more difficult owing to the vast number of cells, cell types, and their interactions.

As a step towards tissue-scale modeling of cells, we consider the problem of simulating large and heterogeneous populations of non-interacting cells. That is, we wish to model the wide variety of cellular states and dynamics experienced by a heterogeneous population of cells in isolation from each other. The output of the approach we propose will be the frequency with which certain events happen over time in a large population. This result could be obtained by averaging a large number of traditional single-cell simulations spanning the entire population, but in extremely heterogeneous populations it becomes infeasible to simulate each distinct subpopulation, in which case the traditional recourse is to estimate the population statistics by Monte Carlo (random sampling) [3]. By design, the basic random sampling procedure captures typical outcomes of these simulations, and only rarely finds atypical occurrences. Yet some biological processes are determined by outliers [4], such as the initiation and development of cancerous cells [5–7] or immune cell clonal selection [8]. If something is known about the circumstances leading to a rare outcome, one can bias Monte Carlo to oversample that outcome and then correct for the biased sampling (a strategy known as importance sampling [9]), but oversampling inevitably introduces sampling biases.

Here we propose an alternative, exact method for simulating heterogeneous populations, which takes advantage of the observation that discrete models have a finite set of possible states. For these models, one can write the instantaneous state of some individual within the population (i.e. a single instance of a Boolean network) using a vector b_(α) that has an entry for every possible state, where we place a 1 in the position corresponding to state α of the individual and a 0 everywhere else. Assuming deterministic dynamics, the time evolution of this individual can then be written as a linear (matrix) operator F _b , so that the time evolution to a different state β occurs via repeated matrix multiplications of the state vector: \(\mathbf {b}_{(\beta)}(t) = F_{b}^{t-t_{0}} \mathbf {b}_{(\alpha)}(t_{0})\) where t₀ and t are integer time steps. This approach is always possible in principle for discrete systems, even when it is too computationally expensive to be feasible in practice.

The usefulness of a linear representation is that the same equations that simulate an individual automatically generalize to simulating arbitrary mixed populations of individuals in different states. The basic idea, illustrated in Fig. 1, is that a population-averaged vector \(\left < \mathbf {b} \right > = \sum _{\alpha } w_{\alpha } \mathbf {b}_{(\alpha)}\) evolves according to the same time evolution operator F _b as does a vector representing an individual, owing to the superposition property of linear systems. We will exploit this fact when we derive a time evolution operator using an algebra tailored for an individual, and then repurpose that operator to simulate mixed populations.

The obvious drawback of the linear method of time evolution is the typically huge number of states in a given system, causing both the state vector b and the time evolution operator F _b to be infeasibly large. For example, Boolean networks are a class of simple models built entirely from ON/OFF variables, yet even these models have an exponential number of states (2 ^N for N Boolean variables in the model). A linear representation of the dynamics is usually only feasible if one looks at a small subspace of the full linear space – i.e. considers only a small subset of the possible set of states of the model. However, this strategy is incompatible with our interest in simulating massively heterogeneous populations: a simulation involving n _s subpopulations necessarily involves at least n _s nonzero entries in both the state vector for the population <b> and in the time evolution matrix operator F _b . Under this proportional scaling with heterogeneity, a linear representation of the system offers no computational improvement over simple enumerative simulation of each individual subpopulation.

Our proposed solution to the problem of proportional scaling with heterogeneity is to change linear basis from the state basis b to a new basis x in which the size of the subspace of interest scales independently of the heterogeneity of the mixed population. (The full spaces of b and x are necessarily of the same size.) Specifically, the goal of this paper is to introduce a linear basis for Boolean variables that we term a ‘product basis’, and give a prescription for calculating the time evolution operator in this basis. Within the Boolean framework, the product basis method is very general, applying to deterministic Boolean models [10], as well as probabilistic [11] and continuous-time [12] Boolean models in the large population limit. Note that although we assume a synchronous updating rule for all variables in the model, this is without loss of generality because an asynchronous network can be modeled as a probabilistic network with synchronous time steps [13].

By way of comparison, we note the significant advances that have been made in analyzing the possible long-term outcomes, or attractors (steady states and limit cycles), of Boolean models [10, 14–22]. Attractors have been found using network reduction algorithms that find simple networks encoding the long-term behavior of more complex networks [14, 16, 23], methods that solve steady states as zeros of a polynomial equation [24], SAT methods [18, 19, 25], and binary decision diagrams [20, 21, 26] (the introduction of [15] provides a review of these techniques). These techniques differ from our proposed method in that they focus on possible long-term behaviors, whereas ours gives an explicit population-averaged simulation of a defined starting population. While our method does have the ability to generate long-term dynamical equations that can be used to find attractors (see Additional file 1: Appendix 2), it differs in that it finds the attractors of a set of variables of interest, not the attractors of the complete state.

We first derive the procedure for calculating a product basis simulation to track transient and long term model behaviors. Initially we consider simulations of individuals, defined throughout this work as single instances of a Boolean network; then we use the same mathematics to directly simulate mixed populations of individuals in many different states. A toy example calculation is given to illustrate the mathematics. Finally, we demonstrate an application of our method by using it to simulate a large heterogeneous population of individuals whose dynamics is described by a published T-cell network [27].

Here we consider a Boolean network model that consists of N variables, which updates its state at each time step using a deterministic update rule. Initially we will focus on individuals described by this model, whose instantaneous state is described by Boolean values for model variables B₁,B₂…B _N which evolve according to the rules of the model. We index these model variables using Roman letters (e.g., i,j,…), and use Greek letters (e.g., α,β,…) to refer to subsets of the model variables. For each possible subset of model variables α={i,j,… }, there is a unique state basis variable b _α (which we sometimes write as b_ij…) and a unique product basis variable x _α =x_ij…, which we formally define below. We consider the state basis variables to be the components of a vector b, and the product basis variables to be components of vector x. There are 2 ^N subsets of all N Boolean model variables (including both the set θ of all model variables and the empty subset denoted ∅), so both b and x are vectors of length 2 ^N .

The empty product x _∅ equals 1. The fact that each x _α is a simple product of model variables motivates the terminology ‘product basis’. As shown in Additional file 1: Appendix 1, a product basis vector x is a fully equivalent representation of a mixed population to the state basis vector b.

Simulations of mixed populations

We build our simulations by selecting a set of product basis variables of interest and associating an update rule f _α with each variable x _α so that x _α (t+1)=f _α (x(t)) (the exception being the case of continuous-time Boolean networks in which f _α =dx _α /dt, but we will treat those separately later). We construct the simulation in two steps. The first step is to build the single-index update rules f _i (i.e. α={i}) over all model variables i, by enumeration of their input states. The second step is to build certain multi-index update rules f_ij… as needed until the system of equations closes (i.e. until we have solved for each f _α corresponding to a variable x _α appearing in another time evolution equation f _β ). To begin with, we show how to build the change-of-basis operator T that converts state space basis vectors to product space basis vectors through the formula x=T·b.

Algorithm 1

(Constructing a change-of-basis matrix) Consider any set θ^′ containing n≤N of the model variables, for which κ and α denote subsets. Define T⁽ⁿ⁾ as the change-of-basis matrix that converts a length- 2 ⁿ b vector indexed by κ to a length- 2 ⁿ x vector indexed by α, and let \(T^{(n)}_{\alpha \kappa }\) denote the matrix element projecting b _κ onto x _α . We construct T⁽ⁿ⁾ by assigning a 1 to each matrix element \(T^{(n)}_{\alpha \kappa }\) for which α⊆κ, and 0 to all other elements.

Proof

Consider the state vector b_(κ) whose entries are all 0 except for a 1 in the position corresponding to state κ, which describes an individual in state κ. The product basis representation x of this individual is found by multiplying \(x_{\alpha } = \sum _{\gamma \subseteq \theta '} T^{(n)}_{\alpha \gamma } b_{\gamma } = T^{(n)}_{\alpha \kappa }\); thus x equals the column of the change-of-basis matrix that multiplies b _κ . Using Definition 1, the value of the product basis variable x _α , corresponding to matrix element \(T^{(n)}_{\alpha \kappa }\), is 1 if α⊆κ and 0 otherwise. □

We can now provide a procedure for calculating the single variable update rules f _i . To do so, we consider only the relatively small subset of variable i’s inputs, rather than the full set of model variables. We use a superscript [ i] to denote quantities pertaining to the input subset; thus we define N^[i] as the number of inputs to model variable i, θ^[i] as the set of those input variables, b^[i]={b _ρ ∣ρ⊆θ^[i]} as the state space of input variables, and x^[i]={x _ρ ∣ρ⊆θ^[i]} as the corresponding product space. In biological networks, N^[i] is usually small enough that we can explicitly write the change-of-basis operator \(T^{\left (N^{[i]}\right)}\) in this space using Algorithm 1.

Algorithm 2

(Computing f _i ) Define k^[i] as a row vector such that \(k_{\alpha }^{[i]}\) is 1 if the pattern of Boolean inputs \(b_{\alpha }^{[i]}\) produces a 1 in output variable i, and 0 otherwise. Then \(f_{i} = \mathbf {k}^{[i]} \left (T^{\left (N^{[i]}\right)} \right)^{-1} \mathbf {x}^{[i]}\), which is a linear equation in x^[i]⊆x.

Proof

By definition f _i =k^[i]·b^[i], as this expression reproduces the output rule. Using the fact that \(T^{\left (N^{[i]}\right)}\) is invertible (proved in Additional file 1: Appendix 1), we write \(f_{i} = \mathbf {k}^{[i]} \left (T^{\left (N^{[i]}\right)} \right)^{-1} T^{\left (N^{[i]}\right)} \mathbf {b}^{[i]}\) and note that \(T^{(N^{[i]})} \mathbf {b}^{[i]} = \mathbf {x}^{[i]}\), which proves the formula. □

Using the set of single-index f _i , one can compute the linear time evolution function of any multi-index product basis variable f_ij… using the following method.

Algorithm 3

(Computing f_ij…) First compute f _α for α={i,j,… } as f _α ←f _i (x)·f _j (x)·… (expressed in terms of x-basis variables). Next, distribute each term inside the product, so that f _α is a weighted sum of products of x-basis variables, and replace each nonlinear product of terms x _β ·x _γ ·… appearing inside f _α with the product basis variable x _μ where μ=β∪γ∪⋯. This gives an expression for f _α that is linear in x.

Proof

First we show that f _α =f_ij… equals the product f _i ·f _j ·…:

$$\begin{array}{*{20}l} f_{ij\dots}(t) &= x_{ij\dots}(t+1)\\ &= x_{i}(t+1) \cdot x_{j}(t+1) \cdot \dots\\ &= f_{i}(t) \cdot f_{j}(t) \cdot \dots\\ &= f_{i} \cdot f_{j} \cdot \ldots \end{array} $$

where the last line emphasizes that there is no time dependence in f_ij….

The second step is to show that each x _β ·x _γ ·… equals x_{β∪γ∪…}. Let k,l,… be the elements of μ=β∪γ∪…; then \(x_{\beta } \cdot x_{\gamma } \cdot \ldots = x_{k}^{p_{k}} \cdot x_{l}^{p_{l}} \cdot \dots \) where p _k ,p _l ,⋯≥1 are the respective number of times k,l,… appear in β,γ,⋯. Since all x _i are Boolean variables, we have \(x_{i}^{p} = x_{i}\) for any p≥1. Thus x _β ·x _γ ·…=x _k ·x _l ·…=x _μ . □

Using Algorithms 2 and 3, we can give the full procedure for building a simulation that time evolves any arbitrary set of product basis variables of interest describing some individual modeled by the Boolean rules. We denote the set of product variables of interest as Ω₀; note that each element of Ω₀ is itself a set of indices over model variables. Our algorithm constructs f _α for each α∈Ω₀, then additional f _β to evolve each x _β that appears in the formula for f _α , etc. until the equations form a closed system (i.e. there is an update rule for every product basis variable appearing in the system of equations).

Algorithm 4

(Building a product basis simulation) First compute the full set of single-index f _i using Algorithm 2. Next, initialize the total set of required product basis variables Ω←Ω₀, and the set of product basis variables whose dynamics we have solved Ω _S ←∅. Finally, iteratively solve f _α for each α∈Ω∖Ω _S using Algorithm 3, while updating Ω _S ←Ω _S ∪{α} and Ω←Ω∪{μ} for each variable x _μ appearing in the formula for f _α . Iteration continues until Ω _S =Ω.

Proof

We first note that Ω eventually converges to a finite set owing to the fact that the number of variables |θ| is finite, and that Ω⊆θ never loses elements of θ at each iterative step. Since a) Ω _S ⊆Ω, b) Ω _S accumulates one term in Ω∖Ω _S at each iterative step, and c) Ω converges, the algorithm always ends in a finite number of steps. When the algorithm terminates, the set of solved variables Ω _S equals the set of variables Ω appearing in the equations, so the resulting system of equations is closed. □

Writing our final closed system of linear equations as a square matrix F, we have a very simple update rule for simulations using the product basis variables: x(t+1)=F·x(t).

Our final step is to generalize to a mixed-population simulation.

Definition 2

Define a population-level state vector <x> as a population-weighted linear combination of the state vectors of the subpopulations x_(α):

$$\begin{array}{*{20}l} \left< \mathbf{x} \right>(t) &:= \sum_{\alpha} w_{\alpha} \mathbf{x}_{(\alpha)}(t) \end{array} $$

(1)

where the weighting factors w _α are proportional to the occurrence of subpopulations α in the overall mixed population.

Note that we are free to choose the normalization constant \(W = \sum _{\alpha } w_{\alpha }\). In all of our examples, we will choose W=1, leading to the interpretation that w _α is the fraction of the population in state α. Irrespective of the choice of W, this representation can be used to evolve mixed populations over time using our existing operators.

Claim 1

Each element <x _α > of the vector <x> is proportional to the overall occurrence of individuals having x _α =1 in the mixed population. For W=1, <x _α > is the fraction of the population having x _α =1.

Proof

Equation 1 is the definition of a weighted average of x, whose weighting factors w are proportional to the population fraction. □

Claim 2

Our linear operator F that evolves any arbitrary state of an individual x _α over time also correctly evolves the state of any arbitrary mixed population <x> over time.

Proof

This Claim follows from the fact that F commutes with the sum over subpopulations:

$$\begin{array}{*{20}l} \left< \mathbf{x} \right>(t) &= \sum_{\alpha} w_{\alpha} F^{t-t_{0}} \mathbf{x}_{(\alpha)}(t_{0})\\ &= F^{t-t_{0}} \sum_{\alpha} w_{\alpha} \mathbf{x}_{(\alpha)}(t_{0})\\ &= F^{t-t_{0}} \left< \mathbf{x} \right>(t_{0}). \end{array} $$

(2)

□

Thus, the time evolution operator F produced by Algorithm 4 correctly evolves the mean occurrence of variables <x> over time in any mixed population, despite the fact that this operator was derived for an x that describes an individual (notably in assuming that each x _α is Boolean).

Example 1: building equations

Consider the model shown in Fig. 2, whose Boolean model variables update according to the following rules:

$$\begin{array}{*{20}l} A(t+1) &= \overline{B(t)}\\ B(t+1) &= A(t) \land C(t)\\ C(t+1) &= A(t) \lor B(t). \end{array} $$

To build a product basis simulation, we first compute the change-of-basis matrices that will be used to compute the single-variable update rules f _A , f _B and f _C . Variable A takes input from the single variable B, so calculating f _A requires the change-of-basis matrix T⁽¹⁾. Ordering the elements by the subscripts (∅,B) respectively, and applying Algorithm 1, we obtain

$$\begin{array}{*{20}l} T^{(1)} &= \left[ \begin{array}{cc} 1 & 1\\ 0 & 1\\ \end{array} \right]. \end{array} $$

Variables B and C each take input from two variables. Ordering B’s input states (∅,A,C,AC), and C’s input states (∅,A,B,AB), we find that both f _B and f _C are computed using

$$\begin{array}{*{20}l} T^{(2)} &= \left[ \begin{array}{cccc} 1 & 1 & 1 & 1\\ 0 & 1 & 0 & 1\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1\\ \end{array} \right]. \end{array} $$

Next, we build the single-index update functions. Variable A takes input only from variable B, so the possible patterns of active inputs are (∅,B)^⊤, corresponding to the state basis variables \(\left (b^{[A]}_{\emptyset }, b^{[A]}_{B}\right)^{\top }\). The respective outputs are (1,0)=k^[A] due to the NOT gate, from which we can immediately calculate f _A using Algorithm 2:

$$\begin{array}{*{20}l} f_{A} &= \mathbf{k}^{[A]} \left(T^{(1)} \right)^{-1} \left[ \begin{array}{c} x_{\emptyset}\\ x_{B} \end{array} \right]\\ &= \left[ \begin{array}{cc} 1 & 0 \end{array} \right] \left[ \begin{array}{cc} 1 & -1\\ 0 & 1 \end{array} \right] \left[ \begin{array}{c} x_{\emptyset}\\ x_{B} \end{array} \right]\\ &= x_{\emptyset} - x_{B}. \end{array} $$

(3)

In the same way we find that the input patterns \(\left (\emptyset, b^{[B]}_{A}, b^{[B]}_{C}, b^{[B]}_{AC}\right)^{\top }\) to variable B lead to outputs k^[B]=(0,0,0,1), and the input patterns \(\left (\emptyset, b^{[C]}_{A}, b^{[C]}_{B}, b^{[C]}_{AB}\right)^{\top }\) to variable C lead to outputs k^[C]=(0,1,1,1). Using these together with T⁽²⁾ we compute:

$$\begin{array}{*{20}l} f_{B} &= \left[ \begin{array}{cccc} 0 & 0 & 0 & 1 \end{array} \right] \left[ \begin{array}{cccc} 1 & -1 & -1 & 1\\ 0 & 1 & 0 & -1\\ 0 & 0 & 1 & -1\\ 0 & 0 & 0 & 1\\ \end{array} \right] \left[ \begin{array}{c} x_{\emptyset}\\ x_{A}\\ x_{C}\\ x_{AC} \end{array} \right]\\ &= x_{AC}\\ f_{C} &= \left[ \begin{array}{cccc} 0 & 1 & 1 & 1 \end{array} \right] \left[ \begin{array}{cccc} 1 & -1 & -1 & 1\\ 0 & 1 & 0 & -1\\ 0 & 0 & 1 & -1\\ 0 & 0 & 0 & 1\\ \end{array} \right] \left[ \begin{array}{c} x_{\emptyset}\\ x_{A}\\ x_{B}\\ x_{AB} \end{array} \right]\\ &= x_{A} + x_{B} - x_{AB}. \end{array} $$

(4)

Having built the single-index update functions, we can now derive a linear system that tracks the time evolution of any set of product variables that we aim to follow. Suppose we wish to follow the time evolution of variable x _A , corresponding to the activity (or mean activity at the population level) of the Boolean model variable B _A . The immediate equation for this purpose is f _A , which we already derived (Eq. 3), but since it involves x _B and x _∅ , our simulation must also track those variables through time using Eq. 4 along with

$$\begin{array}{*{20}l} f_{\emptyset} &= x_{\emptyset}. \end{array} $$

(5)

Equation 4 requires that we track a new multi-index variable x _AC , requiring us to solve f _AC using Algorithm 3:

$$\begin{array}{*{20}l} f_{AC} &= (1-x_{B}) \cdot (x_{A} + x_{B} - x_{AB})\\ &= x_{A} + x_{B} - x_{AB} - x_{AB} - x_{B} + x_{AB}\\ &= x_{A} - x_{AB}. \end{array} $$

(6)

We continue the process of identifying new variables and solving for their update functions until the equations form a closed system:

$$\begin{array}{*{20}l} f_{AB} &= (1-x_{B}) \cdot (x_{AC})\\ &= x_{AC} - x_{ABC} \end{array} $$

(7)

$$\begin{array}{*{20}l} f_{ABC} &= (1-x_{B}) \cdot (x_{AC}) \cdot (x_{A} + x_{B} - x_{AB})\\ &= x_{AC} - x_{ABC}. \end{array} $$

(8)

Equations 3-8, together with initial values for x _A , x _B , x _AC , x _AB and x _ABC , describe the time evolution of these quantities in an individual Boolean network as a sequence of 0s and 1s in each variable. The final step is to reinterpret these equations as describing the dynamics of a mixed population, formally by taking the mean of both sides of each equation:

$$\begin{array}{*{20}l} f_{\left< A \right>} &= \left< x_{\emptyset} \right> - \left< x_{B} \right>\qquad\qquad\qquad\qquad\qquad\quad \text{(3{{b}})}\\ f_{\left< B \right>} &= \left< x_{AC} \right>\qquad\qquad\qquad\qquad\qquad\qquad\quad\;\, \text{(4{b})}\\ f_{\left< \emptyset \right>} &= \left< x_{\emptyset} \right>\qquad\qquad\qquad\qquad\qquad\qquad\qquad \text{(5{b})}\\ f_{\left< AC \right>} &= \left< x_{A} \right> - \left< x_{AB} \right>\qquad\qquad\qquad\qquad\qquad\;\; \text{(6{b})}\\ f_{\left< AB \right>} &= \left< x_{AC} \right> - \left< x_{ABC} \right>\qquad\qquad\qquad\qquad\quad\;\; \text{(7{b})}\\ f_{\left< ABC \right>} &= \left< x_{AC} \right> - \left< x_{ABC} \right>.\qquad\qquad\qquad\qquad\quad\;\, \text({8{b})} \end{array} $$

The angle-brackets denote an average, and we have used the notation <x _i (t+1)>=f_<i>. Per Claim 2, the linear equations are unaffected by the averaging process, so the same equations used to derive the dynamics of an individual also describe the mean values of those same variables in a mixed population. Whereas the product basis variables take on binary values for an individual, the population-averaged variables are real-valued on the interval [ 0,1] (using our recommended normalization in Claim 1). For example, we would set <x _A >=0.4 if gene A is ON in 40% of the population.

Probabilistic and asynchronous Boolean networks

The product basis method can be applied to probabilistic Boolean networks (PBNs) [11, 28], in which several state transitions are possible at each time step with different probabilities. As we will show, our algorithm works in the large-population limit for which time evolution of the average state <x> is essentially deterministic, despite the fact that each individual PBN is stochastic.

Applying our method to PBNs requires that we reinterpret the meaning of the time evolution equations. For an individual we write:

$$\begin{array}{*{20}l} p(\mathbf{x}(t+1)) = F \cdot \mathbf{x}(t) \end{array} $$

(9)

where p(·) denotes the probability of an outcome. The product basis method still works with this new interpretation of the time evolution operator F, although we note several changes to the logic. First, in Algorithm 2 we generalize each k^[i] to be a vector of likelihoods that each respective input pattern produces a 1 in the output variable, so that f _i =k^[i]·b^[i] as before. Second, the multiplication rule used in Algorithm 3 still holds if updated to read \(p_{x_{ij\dots }} = p(x_{i}) p(x_{j}) \dots \), owing to the independence of outcomes p(x _i ), p(x _j ), etc. Third, we point out that although p(x) is real-valued, the state of an individual x is still binary, so \(x_{\alpha }^{p \ge 1} = x_{\alpha }\) as before (Algorithm 3).

Despite the fact that our algorithm produces valid product basis equations of the form of Eq. 9, the resulting linear system of equations does not form a closed system, simply because the left-hand side uses different variables than the right (p(x) versus x). Our method therefore cannot be used to simulate an individual instance of a PBN. However, by averaging both sides and taking the large-population limit so that p(x)→<x>, the system closes and reproduces Eq. 2. Thus, our system of equations accurately tracks the deterministic dynamics of arbitrary mixed populations of PBNs in the infinite-population limit, despite being unable to model stochastic individuals.

Large populations of asynchronous networks behave identically to large populations of PBNs [13] if we define a uniform time step: the probabilities of the various possible updates over that time step give the state transition weights in the corresponding synchronous PBN. If this time step is small enough, then the likelihood of two causally-connected asynchronous updates happening in the same step is small, and in this limit the local update rules for a PBN accurately model the asynchronous network. Therefore our analysis also applies to large populations of asynchronous networks for small time steps.

Continuous-time Boolean networks

Probabilistic networks give one way to incorporate rate information into our model; another way is to work in continuous time using differential equations [12]: f _A =dx _A /dt. The differential form does require one change in our method: the rate of change of a higher-order variable is found by using the product rule of derivatives. Whereas under a discrete update f_ABC… is the product f _A ·f _B ·f _C ·…, for the differential case we compute:

$$\begin{array}{*{20}l} f_{ABC\dots} &= \frac{d}{dt} \left(x_{A} x_{B} x_{C} \dots \right)\\ &= x_{B} x_{C} \dots f_{A} + x_{A} x_{C} \dots f_{B} + \cdots. \end{array} $$

(10)

Also, under discrete updates, the trivial function is f _∅ =1, but with differential updates it is f _∅ =0.