The main steps of the method adopted in this report are as follows:

1) Selection of breast cancer-related genes and subsequent gene regulatory network construction based on this gene set.

This procedure allowed us to highlight attractors and related genes constantly expressed in the dataset of different patients.

This section consists of 6 subsections. In subsection Choice of the Elements of the Gene Regulatory Network, we describe the procedure by which we selected the constituent elements of the gene regulatory network used in this report. The description of Boolean formalism used to model the network dynamics is in subsection Construction of the Boolean Network Model. In subsection Single-Cell RNA-Seq Data we describe the scRNA-seq data used to quantify the network genes. In subsection Binarization of scRNA-Seq Data, we illustrate the method by which the scRNA-seq values assigned to the constituent elements of the gene network have been binarized, and describe the tool used in this report to obtain this result. The last subsection (Search for Attractors) describes the essential characteristics of the network’s attractors and the procedure, through an appropriate software tool, of its identification by simulating trajectory dynamics.

The following subsections detail the steps shown in Figure 1.

Hallmarks of cancer represent groups of acquired biological features that are critical for its development (Hanahan and Weinberg, 2011). We considered two of these hallmarks, UNLIMITED REPLICATIVE POTENTIAL, and EVASION OF CELL DEATH, as starting points for constructing a representative gene regulatory network of cancer. This modeling strategy was chosen to reduce cancer cell proliferation and promote their death. We then obtained four lists of genes from the MSigDB repository (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp) based on the two hallmarks previously considered, each list representing a specific cellular pathway. The gene lists related to Apoptosis and TP53 represent the “Evasion of cell death” hallmark (Wong, 2011), while Kras pathway (up and down-regulated genes) is indicative of the “Unlimited replicative potential” hallmark (Jančík et al., 2010; Aubrey et al., 2016). The choice of these pathways is justified by their particular relevance in the formation and evolution of tumors, along with the potential model scalability. We retained only the significantly differentially expressed genes (DE) between the MDA-MB-231 cell line, a metastatic triple-negative breast cancer subtype (TNBC), and the MCF10A cell line, used as control (Carels et al., 2015).

The selected genes were analyzed considering the number of interactions (edges) of their respective proteins (vertices) in the interactome. The human interactome used in this report is from the intact-micluster.txt file (IntAct database, version updated December 2017 accessed on January 11, 2018, at ftp://ftp.ebi.ac.uk/pub/databases/intact/current/psimitab/intact-micluster.txt). Proteins with edge numbers equal to or greater than 50 were selected as seeds to build the gene regulatory network. Those proteins are potential hubs, for which inhibition has been widely associated with regulatory network disruption (Carels et al., 2015).

We also added five genes to the analysis (HSP90AB1, YWHAB, VIM, CSNK2B, and TK1) (Carels et al., 2015) whose knockdown was shown to inihibit the cell growth and promote the cell death of MDA-MB-231 in vitro (Tilli et al., 2016).

We used the human interactome to define the connections between the proteins coded by the selected genes (network vertices). In case of a lack of a direct relationship between two vertices, we looked for possible intermediary vertices (up to three). We excluded intermediary vertices absent in the gene expression data or with low expression values in MDA-MB-231.

We enriched this network with transcriptional factors that regulate the selected vertices, i.e., differentially expressed hubs and intermediary proteins. We performed this analysis with the online tool TRRUST (Han et al., 2015).

The human interactome from IntAct defines the direction of the interactions (node A regulates node B), but not their function (activation or inhibition). For the definition of interaction functions, we used the Metacore algorithm (Ekins et al., 2007).

We constructed a directed graph model based on Boolean logic from scRNA-seq data. The vertices represent the constituent elements of the dynamic cellular model, and their connections are for the functional regulations acting between them (Emmert-Streib et al., 2014). Boolean network modeling is among the simplest methods for dynamic modeling (Thomas, 1973), but at the same time with characteristics of reliability in providing insights into the dynamics of a system (Herrmann et al., 2012; Siegle et al., 2018). We have translated the gene expression status of a gene into the value of a Boolean variable ( B ) , which can be True or False (1 or 0) based on RNA-seq data. Thus, for the n vertices of our network, we have: X = { x 1 ,   x 2 ,   x 3 , … . , x n }   ,   x i ∈ B (1) This formalization finds its justification considering that one can describe many biological processes, such as concentration levels, through the Hill-Function. For most of the Hill function coefficient values, the resulting curve is a sigmoidal curve, which can be approximated by a dichotomous step-function (Schwab et al., 2020).

The representation of this network’s state in the discrete-state flow of time is a vector whose components are the network’s vertices: x → = ( x 1 ( t ) , … … , x n ( t ) ) (2) and the passage from a certain point of the state space of the system to another is due to the regulatory action of the corresponding Boolean functions: x i ( t + 1 ) = f i ( x → ( t ) ) , f i : B n → B (3) for n nodes of the network.

We decided to adopt a synchronous update mode, where all genes update their values simultaneously at consecutive time points: T ( x → i t ,   x → i t + 1 ) = T 1 ( x 1 t , x 1 t + 1 ) ∧ … . ∧ T n ( x n t , x n t + 1 ) (4) In the Eq. 4, where T ( x → i t ,   x → i t + 1 ) represents the transition function of the state of the network, all the genes in the network simultaneously perform the transition from the state x → i t to the next state     x → i t + 1 in transitions T 1 , T 2 , … , T n occorring in the system. Some reports state that asynchronous updating seems better to model biological systems (Schwab et al., 2020). Nevertheless, synchronous dynamic evolution is computationally more efficient for the type of network used in this report and seems to represent the network’s dynamic behavior in a very similar way (Schwab et al., 2020).

Identifying the rules of interaction among the different entities of a network is usually one of the most challenging tasks in studying gene regulatory network systems. Our choice was oriented towards the nested canalyzing functions (Hinkelmann and Jarrah, 2012), where multiple variables act simultaneously on the function, determining a mechanism of domination of one variable or a group of variables concerning the others based on their Boolean state. For example, in the expression (A ∧ B) ⋁ (C ∧ D), if A ∧ B = 1, the first two variables dominate and determine the expression value. If (A ∧ B) = 0, the expression value is defined by (C ∧ D). Furthermore, it has been shown that nested canalizing functions are a good representation of biological regulations (Nikolajewa et al., 2007) (Harris et al., 2002).

The scRNA-seq data were obtained from the NCBI Gene Expression Omnibus database (accession number GSE 75688, accessed in March 2020). These data refer to the genomic expression profile of 11 patients with 549 cells analyzed. Most of those cells were malignants, while others were not. A large part of the latter were immune T-cells, immune B-cells, and myeloid immune cells. The cancer cells analyzed represented the four subtypes of breast cancer: luminal A, luminal B, HER2, and TNBC (Chung et al., 2017). We used single-cell data for the analyzed network’s corresponding genes, excluding data related to pooled samples (bulk RNA-seq) (Supplementary Table S1). The four subtypes of breast cancer were present among the samples of the 11 patients: BC01_X and BC02_X for luminal A, BC03_X for luminal B, BC04_X, BC05_X and BC06_X for HER2, BC07_X, BC08_X, BC09_X, BC10_X and BC11_X for TNBC. For BC03_X and BC07_X patients, there were metastatic lymph datasets corresponding to BC03lN and BC07LN. For the patient BC09_X, there was another single-cell RNA-seq (BC09Re_X). Note that patient BC05 is the only patient who received prior treatment (neoadjuvant chemotherapy and Herceptin).

As specified in the above description, the different types of breast cancer encountered in this report have already been identified in the dataset.

It is worth pointing out that the model is associated with a specific group of cells for each patient. Nevertheless, this approach can be conceptually made equivalent to the one defined as a multi-cell pathway, and that the relatively high number of available cells analyzed allowed a correct use of the R “Binarize” application, used in this report for the extraction of the Boolean value.

Once the genes making up the network were found and its topology defined, and finally assigned the corresponding scRNA-seq values to each element of the network, the next operation necessary for the Boolean network modeling of the system was to binarize the gene expression values assigned to each single node, such as f : ℝ → B using f ( u ) = { 0 u ≤ t 1 u ≥ t (5) where t is the separation threshold. This result was achieved through the use of the BASC-B algorithm (Binarization Across multiple SCales) (Hopfensitz et al., 2012). The BASC algorithm considers as input values a sorted vector in ascending order ( u 1 , … , u N ) ∈ ℝ , and based on it, BASC defines a discrete, monotonically increasing step function f(x) with N steps and N - 1 discontinuities: f ( x ) = ∑ i = 1 N u i I A i ( x ) (6) with i ∈ { 1 , … , N } . Defining d i = N − 1 as discontinuities, we have A i as intervals defined as follow A i = { ( 0 , d i ] , i f   i = 1 ( d i − 1 , N ] , i f   i = N ( d i − 1 , d i ] , o t h e r w i s e (7) and I A as I A i ( x ) = { 1 , i f   x ∈ A 0 , o t h e r w i s e (8) Once the step function f(x) is obtained using the output vector ordered in increasing order, the algorithm calculates additional step functions that approximate this function with a smaller number of discontinuities. The algorithm then finds the strongest discontinuity in each step function and estimates the strongest discontinuities’ location and variation. This algorithm was implemented through the R Software package BiTrinA (Müssel et al., 2016).

After defining the binarized RNA-seq values on each node of the network and establishing the rules that determine its dynamics, we sought the network’s stable equilibrium state, i.e., the attractors, which can be either singleton (composed of a single state) or cyclic (composed of multiple states) (Huang et al., 2009). The hypothesis under which one may consider the malignant state as a particular type of attractor (Huang et al., 2009; Creixell et al., 2012; Yu and Wang, 2016; Poret and Guziolowski, 2018) has oriented our investigations towards the localization and characterization of attractors in Boolean networks. Furthermore, basins of attraction include all the system states that evolve into a given attractor. They conceptually represent the epigenetic barriers that delimit the basin of attraction (Conforte et al., 2020). We obtained the corresponding attractors matching the gene network for each scRNA-seq dataset of the eleven patients with breast cancer (Chung et al., 2017). Attractor analysis allowed us to highlight the key genes in each basin of attraction and how their inhibition could determine a change in cell fate by using the python Open Source software application “BooleanNet” (Albert et al., 2008).

We performed the following procedure to search for attractors from the available data:

• We used BooleanNet (Albert et al., 2008) to assess the logic functions assigned to each gene of our regulatory network and search for Boolean attractors.

Figure 2 summarizes the adopted procedure.

The process of choosing the elements (genes) constituting the gene regulatory network adopted in this report produced the following results.

First, we obtained 761 genes derived from the Broad Institute repository, divided into four lists related to the two cancer hallmarks used to build the network. The 761 genes were classified as follows (Supplementary Table S2):

• 161 related to the APOPTOSIS pathway

In order to retain only differentially expressed genes, we compared the lists obtained with the RNA-seq data of the MDA-MB-231 and MCF10A cell lines (Tilli et al., 2016), obtaining the following results:

• Because they were neither present in the gene expression data of MDA-MB-231 nor in the MCF10A one 1) 129 genes were excluded from the APOPTOSIS pathway list, leaving 32 genes;

Among the genes retained, we selected only those that were differentially expressed, which resulted in a total of 51 genes (Supplementary Table S3):

• 18 genes of the APOPTOSIS group, 9 Up and 9 Down;

Based on the number of interactions in the interactome, 15 genes of the 51 were selected, from which 5 (HSP90AB1, YWHAB, VIM, CSNK2B, TK1), considered more relevant for the present research, have been added to the network (Carels et al., 2015). As outlined above, 1) the vertice vertex connections obtained by comparison with IntAct human interactome, 2) the inclusion of intermediate vertices, 3) the enrichment of the network with transcriptional factors that regulate the selected vertices with the online tool TRRUST (Han et al., 2015), and 4) the activation or inhibition of vertex inputs obtained with the Metacore algorithm (Karin, 2006) (Supplementary Table S4), allowed us to obtain a gene regulatory network consisting of 103 vertices (see Figure 3), and whose dynamics were regulated by nested canalyzing functions (Hinkelmann and Jarrah, 2012) (Supplementary Table S5).

The 14 groups of scRNA-seq binarization values from the 11 patients belong to 4 types of breast cancer. They were divided and organized according to the following criterion: 26 single-cell datasets for the patient BC01_X, 56 for BC02_X, 37 and 55 for BC03_X and BC03LN, 59 for BC04_X, 77 for BC05_X, 25 for BC06_X, 51 and 53 for BC07_X and BC07LN, 23 for BC08_X, 29 and 31 for BC09_X and BC09Re, 16 for BC10_X, 11 for BC11_X (see Figure 4). It is worth noting that the values relating to pooled single-cell present in each patient group were excluded from the count.

The gene expression values of every single cell of each patient were matched to the corresponding 103 genes making up the gene regulatory network and subsequently binarized using the BASC-B algorithm (Hopfensitz et al., 2012) (Supplementary Table S6).

For every single cell of the 14 groups representing the 11 patients of breast cancer, the 103 binarized values at each node of the gene regulatory network were processed by the previously described Boolean attractor search procedure (Albert et al., 2008). The attractors obtained for malignant cells, stromal cells, immune B and T-cells, and myeloid cells are thus summarized as follow: 1) BC01_X: 19 malignant attractors, 2 stromal cell attractors, 5 no result; 2) BC02_X: 49 malignant attractors, 7 no result; 3) BC03_X: 15 malignant attractors, 7 immune B-cell attractors, 5 immune T-cell attractors, 10 no results; 4) BC03LN_X: 6 malignant attractors, 35 immune B-cell attractors, 3 immune T-cell attractors, 11 no results; 5) BC04_X: 42 malignant attractors, 3 immune T-cell attractors, 2 immune Myeloid attractors, 12 no results; 6) BC05_X: 74 malignant attractors, 3 no results; 7) BC06_X: 6 malignant attractors, 2 stromal cell attractors, 6 immune B-cell attractors, 11 no results 8) BC07_X: 24 malignant attractors, 4 stromal cell attractors, 3 immune B-cell attractors, 4 immune T-cell attractors, 8 immune myeloid attractors, 8 no results; 9) BC07LN_X: 24 malignant attractors, 19 immune B-cell attractors, 10 no results; 10) BC08_X: 15 malignant attractors, 6 stromal cell attractors, 2 no results; 11) BC09_X: 2 stromal cell attractors, 1 immune B-cell attractors, 7 immune T-cell attractors, 15 immune myeloid attractors, 4 no results; 12) BC09Re_X: 2 stromal cell attractors, 1 immune B-cell attractors, 20 immune T-cell attractors, 6 immune myeloid attractors, 2 no results; 13) BC10_X: 11 malignant attractors, 2 stromal cell attractors, 2 immune myeloid attractors, 1 no results; and 14) BC11_X: 10 malignant attractors, 1 no results.

Based on these results, we decided to exclude patient data BC09_X and BC09Re due to the lack of specific tumor cell attractors (Supplementary Table S7).

The pie-charts in Figure 5 shows these results expressed as a percentage of the total Single-cell RNA-seq datasets available for each patient, highlighting the success rate in the search for specific attractors on tumor cells.

We selected network genes based on tumor cells’ attractors according to the following criterion: each patient’s attractors kept their level of gene expression (or non-expression) constant for a particular gene, formalized respectively with the symbol “True” or “False”. This criterion allowed us to formulate the following considerations on the results obtained:

• BAX is expressed in the attractors of all patients except BC03_X.

These results are summarized in Figure 6.

Interestingly, PLAT is never expressed in patients with breast cancer classified as TNBC. Unlike for patients of Luminal A and Luminal B in which the inactivity of PLAT affects 50% of patients, this characteristic covers all TNBC group cases (Figure 5).

Further considerations concern the comparison between the attractors of malignant cells with other types of cells from the same patient. For example, in the BC07LN_ patient sample, it is interesting to compare EEF1G in malignant and immune B-cells. In the first case, the gene is expressed in only 4.2% of the attractors detected (1/24), while in the second case, it is expressed in 36.9% of the attractors detected (7/19). Considering DDR1, the attractor rate of expression was 60% (9/15) in the patient sample BC08_X. For stromal cells, the attractor level of expression for the same gene is 100% (6/6).