Simulating heterogeneous populations using Boolean models

Background Certain biological processes, such as the development of cancer and immune activation, can be controlled by rare cellular events that are difficult to capture computationally through simulations of individual cells. Information about such rare events can be gleaned from an attractor analysis, for which a variety of methods exist (in particular for Boolean models). However, explicitly simulating a defined mixed population of cells in a way that tracks even the rarest subpopulations remains an open challenge. Results Here we show that when cellular states are described using a Boolean network model, one can exactly simulate the dynamics of non-interacting, highly heterogeneous populations directly, without having to model the various subpopulations. This strategy captures even the rarest outcomes of the model with no sampling error. Our method can incorporate heterogeneity in both cell state and, by augmenting the model, the underlying rules of the network as well (e.g., introducing loss-of-function genetic alterations). We demonstrate our method by using it to simulate a heterogeneous population of Boolean networks modeling the T-cell receptor, spanning ∼ 1020 distinct cellular states and mutational profiles. Conclusions We have developed a method for using Boolean models to perform a population-level simulation, in which the population consists of non-interacting individuals existing in different states. This approach can be used even when there are far too many distinct subpopulations to model individually. Electronic supplementary material The online version of this article (10.1186/s12918-018-0591-9) contains supplementary material, which is available to authorized users.

Proof. Our proof first requires us to explicitly write the change-of-basis matrix T . We can do this using Algorithm 1 once we have chosen a particular ordering of the entries of both the b and x vectors. Here we will choose to order the entries of both b and x by binary order of the subscripts on their respective state variables. In other words, the index is given by the binary number having 1s in the positions indicated by the subscripts of each state variable: the state variable written as b α = b ij... is the mth entry in the state space vector b and x α = x ij... is the mth entry in the product space vector x, where m = 2 i−1 + 2 j−1 + · · · + 1, for 1 ≤ {i, j, . . . } ≤ N . This choice of ordering fixes a particular representation of the change-of-basis matrix. For example, for the case of three Boolean variables A, B and C, the change-of-basis matrix is For notational convenience, the two lemmas below denote matrix elements using scalar indices as in T mp , rather than using their corresponding set variables as in T αβ . Likewise a product variable x α is sometimes denoted here as x m , referring to the position of x α in x using the above representation.
Our objective is to show that the change-of-basis matrix is always upper-triangular, with 1s on the diagonal.
Lemma 1. The m th diagonal entry T mm equals 1.
Proof. T mm projects b α onto x α for some α, so it can be written T αα . Referring to algorithm 1, α ⊆ α so this matrix element is 1.
Lemma 2. Every entry below the diagonal (i.e. T mp for m > p) equals 0.
Proof. Matrix element T mp projects state basis variable b p = b β onto product basis variable x m = x α . By Algorithm 1 this matrix element is either 0 or 1. Assume it is a 1, implying that α ⊆ β. Then m = i∈α 2 i−1 + 1 and p = i∈β 2 i−1 + 1 = m + i∈(β\α) 2 i−1 , implying that p ≥ m. This contradicts the assumption that m > p, so T mp = 0.
Thus each diagonal entry of the change-of-basis matrix is 1, and the entire lower triangle is 0. Since the change-ofbasis matrix is triangular, the eigenvalues are the diagonal entries, which are all 1. Therefore the determinant is 1, so the change-of-basis transformation is non-degenerate (and volume-preserving) and therefore invertible.

Appendix 2: equation reduction for late time behavior
A mixed-population simulation may or may not be practical, depending on whether the system of linear equations closes with a manageable number of variables. One way to simplify the equations is to focus on the late-time behavior of the system, after quickly-decaying eigenmodes have dropped near zero and time evolution effectively occurs in the smaller subspace of slowly-decaying or persistent eigenmodes. In the limit where we focus on the infinitetime behavior, this subspace is simply the attractors of the system (steady states and/or limit cycles). Our general strategy will be to identify fast-decaying modes as they appear in the time evolution matrix operator F , set them to zero and use the resulting expression as a constraint that can be used to remove variables from later steps of the calculation. In order to be useful for simplifying calculations, these steps must be taken before the linear system closes -in other words, when our partiallycomputed update function F partial is still a non-square matrix.
Our strategy for finding decaying modes in F partial relies on a QR decomposition, which can be performed on non-square matrices. In the case of a deterministic network, all eigenvalues of the full time evolution operator F have modulus 0 or 1 implying that all decaying modes are strict linear dependencies. (Proof : The product space is finite, so every eigenmode λ must either decay to zero in finite time implying that λ = 0, or else repeat in ∆t time steps so that |λ| = 1 and arg(λ) = 2πk/∆t for integer k.) The set of linearly-dependent rows in F is found by identifying (numerically) zero entries along the diagonal of the R matrix from the QR decomposition. The dependencies themselves are then calculated by solving the matrix equation N = KD for polynomial coefficients K, where N is the submatrix that contains the independent rows of F and D is the submatrix that contains the dependent rows of F . For PBNs, the complication is that decaying eigenvalues of the full system are generally nonzero, so the decaying modes are not simple linear dependencies. In order to find modes in the m × n matrix F partial where n ≥ m, we examine the left-hand m × m portion of this matrix which we denote F . Each decaying eigenvalue of the full operator F that can be identified at this stage in the calculation is part of the spectrum of F , although the reverse is not true. We therefore use the spectrum of F as a set of candidate eigenvalues, and for each sufficiently-small candidate λ (based on a userdefined threshold) find the linearly dependent subspace using a QR decomposition of F partial − λI partial , where I partial is another m × n matrix having 1s on the diagonal and 0s elsewhere.
Each constraint is the result of a linear relationship between the update rules appearing after a certain number of time steps, implying that the same linear relationship in the respective product basis variables appears one time step later. The general form of a linear dependency D beginning at time step t c is: where the subscript under the equals sign indicates the range of simulation time over which linear dependency D is valid. There are several ways in which constraints derived from this dependency could help simplify the calculation, but we must be careful that the application of a constraint should not undo the work of a prior constraint, in order to guarantee that the algorithm eventually terminates. Furthermore, the contraint algorithm should have a self-consistent way of removing a variable with two or more simultaneous constraints. Here we give two constraint algorithms that help to simplify late time analyses in the cases of probabilistic and deterministic Boolean models respectively.
Algorithm 5 (Equation reduction method 1). Given linear dependency D, we choose the lowest-index variable x φ (which may be x ∅ = 1) based on a binary representation of its indices φ as in Appendix 1. Next, we remove that variable from all current and future time evolution equations in F by simple substitution: Constraining a time evolution equation valid from time t e using this formula yields a different time evolution equation valid for t ≥ max(t e , t c ).
Proof. The substitution is justified for t ≥ t c , so the resulting time evolution equation is accurate after time t c . A given constraint either removes a variable (if the constraint polynomial contains only one term), or else replaces a variable with higher-index variables. Each replacement with higher-index variables increases the index number of each constrained term by at least 1, and since the maximum index number is 2 N , the process of applying constraints to a given term in the equations eventually completes. Each constraint application permanently removes one variable from the finite 2 N -sized variable space, so the process of constraining the equations also completes, which is only possible if there are no more fastdecaying modes. Thus this process eventually removes the fast-decaying eigenspace, while preserving the slowlydecaying/persistent eigenspace.
We use equation reduction method 1 for PBNs and asynchronous networks as equation reduction method 2 (see Algorithm 6 below) does not apply to these networks. Equation reduction method 1 is relatively inefficient since it removes only one variable per constraint. There is also some indeterminism: if two constraints apply to a given term in an equation, we can choose only one of the two constraints to apply.
Here we demonstrate equation reduction method 1 by using it to find the attractors of the deterministic network from Example 1 in the main text.

Extension of Example 1
After solving Eqs. 1.1-1.6 in Example 1, we find that they contain a single linear dependency: f AB − f ABC = 0. Since the network is deterministic, this dependency guarantees that after the first time step (i.e. t ≥ 2), x AB will equal x ABC , so we write x AB = t≥2 x ABC . We use this fact to remove x AB , giving a new set of steady state equations: Our new set of equations has the same non-zero eigenspace as the original set (1.1-1.6). However, the new equations lack the null eigenspace because we removed the only linear dependency. The 5 × 5 linear system has 5 eigenvalues, all having phases that are multiples of 2π/5, indicating that the sole attractor is a limit cycle with a period of 5 time steps. We acknowledge that in this case the constraint helps find the attractors, but does not simplify the calculation, since it only appears after the full time evolution operator F has been calculated.
Algorithm 6 (Equation reduction method 2). Given linear dependency D, we obtain a constraint from each variable x φ whose indices φ have the property γ\φ = ∅ for every other term x γ in the dependency. Constraining a time evolution equation of a deterministic model valid from time t e involves constraining every variable x ρ in that equation that is factorized by one of the x φ , in the sense that φ\ρ = ∅. Application of the constraint takes x ρ → −(c α /c φ )x α∪ρ −(c β /c φ )x β∪ρ −· · · , yielding an equation that is valid for t ≥ max(t e , t c ).
Proof. For t ≥ max(t e , t c ), both the original time evolution equation and the linear dependency equation are valid, implying that For deterministic models, x ρ = x φ · x ρ since φ ⊆ ρ and indices can be freely duplicated in deterministic models. Substituting the two formulas gives the constraint rule. Since each constraint application either removes x φ or replaces it with variables containing more indices, we can use the same arguments given for equation reduction method 1 to show that the method eventually removes the fast-decaying eigenspace while leaving correct time evolution equations for the slower-decaying modes.
Equation reduction method 2 is considerably more aggressive than equation reduction method 1, because each constraint removes not only one or more variables involved in the constraint, but also any variables that are factorized by those constraints. As such, a dependency involving a variable with a lower index can exclude a significant fraction of the variable space. We demonstrate equation reduction method 2 using Example 2. Example 2: applying constraints using equation reduction method 2 Suppose we want to find the late time behavior of Boolean variables A, C and E in the network from Figure 4. After two iterations of solving for time evolution functions we have: At this point there are two linear dependencies: f AB = f D and f BC = f B , implying that: Since we are only interested in the late time behavior, we will use constraints C1-C3 to simplify the equations. Note that if we were to retain, for example, x ABC , it would be affected by the constraints involving variables x AB , x B , and x BC , thus it would not be possible to enforce all of the constraints by substitution because there is only one B-index on x ABC to substitute for. However, we can simultaneously enforce all constraints by multiplying them: Viewed another way, constraint C1 attaches indices of AB to every variable containing a D, and an index of D to each variable with indices of AB. Constraint C2 adds an index of C to every variable with an index of B.
After constraining our system and removing the time evolution equations of the removed variables, we have: Our new equations require us to solve for another variable (while applying the constraints): x ABCD . (2.8) The system is now closed (5 equations involving 5 variables), so if our goal is to produce a simulation (valid for t ≥ 2) then we are done. However, if our objective is to find the attractors, then we must remove another dependency that was unmasked by the last equation: implies three useful constraints: Each constraint reduces the size of the variable space. For example, constraint C4 removes all variables containing index A with no C, E, or BC. We did not bother to solve for x BC in terms of x A , x B and x C because that gives us a useless constraint: multiplying a variable having BC indices by x C does not change that variable.
After applying constraints C4-C6, we obtain: This produces new dependencies: After applying constraints C7-C9, our example ends with the following closed system of dynamical equations: The attractor is always reached at or before time step 3. The constraint equations (C1-C9) map our variables of interest (x A , x C , and x E ) to linear combinations of variables in the final system. For example, x A is mapped by constraint C4 to x AC + x AE − x ABC , then mapped using constraints C1, C7, and C8 to x AC , whose dynamics are given by the final time evolution equations. The eigenvalues of this final system are -1, 1 and 1, implying a steady state along with a limit cycle having a period of 2 time steps.
When using equation reduction method 2, the order in which we calculate the time evolution of new variables and remove dependencies greatly influences the complexity of the calculation. One general rule is that constraints on the lowest-index variables tend to be most helpful for deterministic networks, because they factorize the largest part of the product variable space. Following this rule of thumb, our implementation solves the time evolution equations for only those variables having the fewest indices between each search for new dependencies. Additionally, our code solves for low-index factors of new variables even if they are not directly involved in earlier equations in hopes that they will produce helpful constraints. In most tests, this prioritization method speeds up the calculations.
Equation reduction method 2 is susceptible to numerical problems that fall into two categories: 1) polynomial coefficients can become very large, and 2) gradual erosion of numerical precision can blur the distinction between zero and nonzero terms. Our code monitors both sources of error; in the latter case, by tracking the floating-point error in each step of the calculation. When either type of calculation error exceeds defined bounds, our method exits with an error message.
The reason polynomial coefficients can become large is that multiplication by constraints can produce several identical terms that add together. For example, the constraint x 1 = x 2 + x 3 takes x 123 → 2x 123 . Typically the variable being multiplied would be found to be zero later in the calculation (so having a coefficient of 2 is not an error in the math), but if many such constraints are applied and the coefficient grows exponentially, then the numerics will eventually fail. We have included a filter that can usually identify these outlier coefficients and remove them using the rule that all individual (not populationaveraged) product variables are either 0 or 1 at all time steps. For example, if f 4 = 2x 123 , we can reason that if f 4 is either 0 or 1, then x 123 cannot be 1 and therefore we have a new constraint x 123 = 0, which simplifies f 4 → 0. In general, we identify these removable terms using the heuristic that if a large coefficient multiplying a combination of integer-weighted variables is greater than the sum of all other coefficients plus one, then that combination of variables (which is an integer for a homogeneous population) must be zero.
Our equation reduction method 2 was tested by running 10 4 experiments comparing Monte Carlo to our product basis method. Each experiment used a unique randomly-generated 10-node Boolean model having 1-4 inputs per node as before, and involved a population of 100 individuals starting in different random states. The comparison was between a) the dynamics of the product basis equations using a random amount of equation reduction (which determines the earliest time point at which the product basis simulation is valid), and b) the population average of the individual simulations. In each case the product basis simulation exactly reproduced the average of the Monte Carlo simulations. These networks were small, but still provided a stringent test of the numerics, as they could theoretically involve linear systems of equations up to size (2 10 ) 2 . During this evaluation, our program also generated additional networks whose numerical error in the polynomial coefficients overstepped a numerical threshold of 10 −6 and were thus omitted from the Monte Carlo comparison; overall 18.5% of tests generated were discarded.