NF-κB fingerprinting reveals heterogeneous NF-κB composition in diffuse large B-cell lymphoma

Introduction Improving treatments for Diffuse Large B-Cell Lymphoma (DLBCL) is challenged by the vast heterogeneity of the disease. Nuclear factor-κB (NF-κB) is frequently aberrantly activated in DLBCL. Transcriptionally active NF-κB is a dimer containing either RelA, RelB or cRel, but the variability in the composition of NF-κB between and within DLBCL cell populations is not known. Results Here we describe a new flow cytometry-based analysis technique termed “NF-κB fingerprinting” and demonstrate its applicability to DLBCL cell lines, DLBCL core-needle biopsy samples, and healthy donor blood samples. We find each of these cell populations has a unique NF-κB fingerprint and that widely used cell-of-origin classifications are inadequate to capture NF-κB heterogeneity in DLBCL. Computational modeling predicts that RelA is a key determinant of response to microenvironmental stimuli, and we experimentally identify substantial variability in RelA between and within ABC-DLBCL cell lines. We find that when we incorporate NF-κB fingerprints and mutational information into computational models we can predict how heterogeneous DLBCL cell populations respond to microenvironmental stimuli, and we validate these predictions experimentally. Discussion Our results show that the composition of NF-κB is highly heterogeneous in DLBCL and predictive of how DLBCL cells will respond to microenvironmental stimuli. We find that commonly occurring mutations in the NF-κB signaling pathway reduce DLBCL’s response to microenvironmental stimuli. NF-κB fingerprinting is a widely applicable analysis technique to quantify NF-κB heterogeneity in B cell malignancies that reveals functionally significant differences in NF-κB composition within and between cell populations.

. NF-κB activation state cannot be determined from NF-κB fingerprints. Computationally simulated NF-κB fingerprints in six cell population specific computational simulations informed by experimental NF-κB fingerprinting ( Figure  4B). Simulation are the same here as in Figure 4C, except for the U2932 cell line where basal IKK activity was increased 100x and the expression of RelA and RelB increased to recapitulate experiment fingerprint ( Figure 4B). 1,000 cells were simulated in each cell population (6,000 simulations in total), with cell-to-cell variability incorporated as described previously (2), cell density is indicated with a contour plot and each cell population is shown in distinct colors. Figure S5. Computational modeling of DLBCL, including receptor-proximal signaling. A) Schematic of the computational model constructed by combining existing models of TLR signaling (3), BCR signaling (4), and NF-κB/IκB regulation (5). All models are run as published, with active IKK species summed from the BCR and TLR models to determine the active IKK input curve to the NF-κB model. Mutations present in DLBCL are indicated in purple. Detail of the combinatorial complexity of NF-κB dimer and inhibitor interactions are omitted.

Supplemental Modeling Methodology
All modeling files are available at https://github.com/SiFTW/NFkBModel. This repository includes the Jupyter notebooks for running the models and producing the figures.
Software and versions Computational simulations were performed using the free DifferentialEquations.jl package, in the free Julia programming language (version 1.6.2). The version number for each package used in this study is provided below: Generating models from reactions, parameters and rate laws. Computational models were encoded across three .csv files: reactions.csv, parameters.csv, and rateLaws.csv (available on GitHub https://github.com/SiFTW/NFkBModel). All parameters in all simulations were kept consistent with published parameters (5), with the exception of cell line-specific parameterizations defined below. Bespoke Python 3 code was used to assemble these files into a single Julia (.jl) file of equations written into a function compatible with the DifferentialEquations.jl package. This python code is available on GitHub (https://github.com/SiFTW/CSV2JuliaDiffEq). Two different collections of these CSV files were used in this manuscript. One defining the NF-κB, and one defining the model in which TLR, BCR and NF-κB signaling are combined. To assemble each of the two models, the CSV2Julia (https://github.com/SiFTW/CSV2JuliaDiffEq) command was run with the appropriate model definition CSV files as arguments to the CSV2Julia function.
Solving models Simulation length, initial conditions, and algorithm used in generating solutions (such as tolerance) are provided in the Jupyter Notebooks https://github.com/SiFTW/NFkBModel). The steady state of each simulation was obtained through an extended simulation from the initial conditions (100,000 minutes), and no species was found to be substantially changing by this time point. The endpoint of this steady state simulation was used as the initial conditions for the time course phase. For conditions where stimulation was required, this was added to the model at the start of the time course (t=0).
Cell-to-cell variability Cell-to-cell variability was approximated (assuming pre-existing differences in expression and degradation rates) as described previously (2). Table S1 shows the parameters which are distributed within the NF-κB model. To produce the simulations in Figure 1 C-D, the rates of transcription for RelA, cRel, and RelB were multiplied by 10, producing the increased RelA, increased cRel, and increased RelB models respectively. Each of these simulations was simulated to steady state followed by a time course response to increased canonical stimulation beginning at t=0. Other than the canonical stimulation, simulated as an increase in NEMO-IKK, all parameters were kept consistent from the time course-phase in the steady state-phase.

Models from gene expression data
In order to produce cell line-specific models from gene expression data, parameters for the rates of transcription of RelA, cRel, and RelB, p50, p100 were adjusted. U2932 and RIVA-specific models were created by adjusting these parameters. Gene expression data was downloaded from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103934) (6). The gene expression (in log10 counts per million reads) was standardized per gene, such that the expression of a particular gene across all cell lines has 0 mean and 1 standard deviation (a z-distribution). The default parameter for expression of each gene was multiplied by 10 z-score . As such a gene with average expression would not have its expression scaled (10 0 = 1), while a gene with expression 1 standard deviation higher than the average would have 10 fold higher expression (10 1 =10). The result was two models, adjusted from the published B-cell model, incorporating cell linespecific gene expression.

Models recapitulating NF-κB fingerprints
To fit the NF-κB model to flow cytometry data from the NF-κB fingerprints RelA and RelB transcription rates were manually adjusted ( Figure 4). In addition to the models in which only RelA and RelB was adjusted, cell line-specific models were created with elevated basal NEMO-IKK activity by multiplying the level of activity during the steady state phase 100 fold. As this resulted in changes to the levels of RelA and RelB, expression of these parameters was adjusted to compensate for the increase in NEMO activity and recapitulate NF-κB fingerprints ( Figure S4).
Creating a model combining NF-κB, TLR and BCR signaling CSV files encoding the reactions, rate laws, and parameters for the NF-κB model (above), were combined with files encoding published TLR and BCR models (3,4). An additional linking module was defined as linking reactions, linking parameters, and linking rateLaws (see the moduleDefinitionFiles