Lower connectivity of tumor coexpression networks is not specific to cancer.

Global level network analysis of molecular links is necessary for systems level view of complex diseases like cancer. Using genome-wide expression datasets, we constructed and compared gene co-expression based specific networks of pre-cancerous tumors (adenoma) and cancerous tumors (carcinoma) with paired normal networks to assess for any possible changes in network connectivity. Previously, loss of connectivity was reported as a characteristics of cancer samples. Here, we observed that pre-cancerous conditions also had significantly less connections than paired normal samples. We observed a loss of connectivity trend for colorectal adenoma, aldosterone producing adenoma and uterine leiomyoma. We also showed that the loss of connectivity trend is not specific to positive or negative correlation based networks. Differential hub genes, which were the most highly differentially less connected genes in tumor, were mostly different between different datasets. No common gene list could be defined which underlies the lower connectivity of tumor specific networks. Connectivity of colorectal cancer methylation targets was different from other genes. Extracellular space related terms were enriched in negative correlation based differential hubs and common methylation targets of colorectal carcinoma. Our results indicate a systems level change of lower connectivity as cells transform to not only cancer but also pre-cancerous conditions. This systems level behavior could not be attributed to a group of genes.

Introduction

Various studies have associated certain network level (topological) features to cancer [1-6]. Studies using protein-protein interaction networks showed that differentially expressed or mutated gene products in cancer are highly connected [1, 2]. Since, protein-protein interaction networks are less condition-specific than gene expression-based networks, these studies give only a stationary view of cancer related genes. Gene co-expression based large scale condition specific networks could be topologically analyzed and compared to reveal possible system-level changes [5, 7]. It was shown that the co-expression networks of colorectal, non-small cell lung, gastric, pancreatic and cervical cancers have lower connectivity as compared to normal samples [5]. Having the same observation for the various cancer types suggests the system level loss of connectivity could be a specific network characteristic of cancer. However, it is unclear how the connectivity changes progress during cancer development. Thus, the presence of datasets for a pre-cancerous stage (adenoma) and a cancerous stage (carcinoma) provides an opportunity for comparing network characteristics of the different phases during colorectal cancer progression at a systems level. Thus, one might gain some insight into the evolution of global gene relationships during the advancement of a non-malignant tumor (adenoma) to a malignant (carcinoma) tumor. Furthermore, to obtain more reliable results in comparative studies, it is important to use matched tissue samples, such that both tumor and adjacent normal tissue samples are obtained from the same individual. We performed a comparative topological analysis based on gene expression values from matched normal-tumor samples.

Results and discussion

Loss of connectivity in colorectal tumor networks

We constructed specific co-expression networks by using gene expression datasets from colorectal adenoma and colorectal carcinoma patient samples. The normal and tumor samples of the datasets were derived from the same patient and were paired. The networks consisted of pairs of highly expressed and varied genes with a significant differential Pearson Correlation Coefficient (PCC) value (see Methods). Several preprocessing steps eliminated genes or gene pairs, such as the genes with low expression levels based on an arbitrary threshold value. The preprocessing steps (differential correlation effect, correlation method; Pearson or Spearman, variation method; standard deviation or quartile coefficient of dispersion) and the threshold values (PCC value threshold, low expression threshold, low variation threshold) had no significant influence on our observations (see Appendix).

We compared the specific networks of the tumors to their respective paired normal networks. Positive and negative gene-gene correlations were analyzed separately. To gain an understanding of the connectivity changes in colorectal cancer, degree values of genes in tumor networks were compared to paired normal networks. The degree values in tumor networks were lower than their paired normal datasets, for both adenoma and carcinoma (Fig. 1). The difference between the normal and tumor was statistically significant (Mann-Whitney paired one-sided test p-value <  10- 145). The same result was observed for both positive and negative correlation based networks. In order to confirm this observation, we analyzed a second independent dataset (see Methods). We observed the same results for the second dataset (Fig. 2, Mann-Whitney paired one-sided test p-value <  10- 200). In addition to degree values, other topological parameters of the tumor and normal networks were also analyzed (Table 1). In general, tumor networks had parameters suggesting lower connectivity. For instance tumor networks had higher characteristic path length (average of shortest path lengths) than normal networks.

Fig.1

Comparison of degree values for Normal (N) and Tumor (T) samples, in colorectal adenoma (Ade 1; GSE8671 dataset) and colorectal carcinoma (Car 1; GSE18105 dataset). Positive (Pos) and Negative (Neg) correlations were analyzed and showed separately.

Fig.2

Comparison of degree values for Normal (N) and Tumor (T) samples, in colorectal adenoma (Ade 2; GSE89076 dataset) and carcinoma (Car 2; GSE89076 dataset). Positive (Pos) and Negative (Neg) correlations were analyzed and showed separately.

Table 1

Statistical topological parameters of Positive Correlation based Normal Networks (PCNN), Positive Correlation based Tumor Networks (PCTN), Negative Correlation based Normal Networks (NCNN), and Negative Correlation based Tumor Networks (NCTN) for colorectal adenoma 1 (Ade 1), colorectal adenoma 2 (Ade 2), colorectal carcinoma 1 (Car 1) and colorectal carcinoma 2 (Car 2) datasets

	Number of Nodes	Number of Edges	Density	Average Degree	Clustering Coefficient	Diameter	Characteristic path length	Number of connected components	Power law value	Power law fit R squared
Ade 1 PCNN	2820	45677	0.01	32.40	0.02	4	2.64	1	0.97	0.26
Ade 1 PCTN	2812	14688	0.004	10.45	0.01	7	3.56	1	1.75	0.65
Ade 1 NCNN	2815	16559	0.004	11.77	0.01	7	3.38	1	1.50	0.65
Ade 1 NCTN	2669	5509	0.002	4.13	0.002	11	5.11	6	2.12	0.86
Ade 2 PCNN	4724	60392	0.005	25.57	0.01	6	2.90	1	1.36	0.56
Ade 2 PCTN	3712	5675	0.001	3.06	0.005	15	5.43	68	1.68	0.79
Ade 2 NCNN	4721	73217	0.01	31.02	0.01	5	2.79	1	1.24	0.46
Ade 2 NCTN	4105	8093	0.001	3.94	0.005	13	5.00	30	1.83	0.83
Car 1 PCNN	7880	863446	0.03	219.15	0.04	3	1.99	1	1.07	0.25
Car 1 PCTN	7880	317777	0.01	80.65	0.02	4	2.39	1	1.48	0.48
Car 1 NCNN	7880	776269	0.03	197.02	0.03	3	2.00	1	1.18	0.31
Car 1 NCTN	7880	222308	0.01	56.42	0.01	4	2.61	1	1.58	0.57
Car 2 PCNN	6090	202178	0.01	66.40	0.02	4	2.44	1	1.25	0.38
Car 2 PCTN	6089	70276	0.004	23.08	0.01	6	3.02	1	1.50	0.59
Car 2 NCNN	6090	139454	0.01	45.80	0.01	4	2.66	1	1.25	0.40
Car 2 NCTN	6089	66193	0.004	21.74	0.01	6	3.08	1	1.52	0.60

We also analyzed the degree distributions of two different pre-cancerous tumors, which were aldosterone producing adenoma (APA) and uterine leiomyoma (UL). Specific coexpression networks for APA and UL were constructed similarly and the degree values of genes in tumor networks were compared to paired normal networks. For both APA and UL, the degree values in tumor networks were lower than their paired normal datasets (Fig. 3). The difference between the normal and tumor was statistically significant (Mann-Whitney paired one-sided tests; APA p-value <  10- 268, UL p-value <  10- 56).

Fig.3

Comparison of degree values for Normal (N) and Tumor (T) samples, in aldosterone producing adenoma (APA; GSE60042 dataset) and uterine leiomyoma (UL; GSE31699 dataset). Positive (Pos) and Negative (Neg) correlations were analyzed and showed separately.

Our results confirmed a previous report by Anglani et al. [5], where the colorectal carcinoma tumor network was significantly less connected than its paired normal specific network. Additionally, we observed the loss of connectivity trend in not only colorectal carcinoma but also in colorectal adenoma, which is considered a pre-cancerous disease stage [8]. We also observed the same trend in two different pre-cancerous tumor types. Our observations suggest the loss of connectivity as a common systems-level trend, and it could be be initiated at an earlier stage of carcinogenesis.

Differentially less connected genes and differential hubs

We observed the same topological change in different datasets which were generated by different microarray platforms (see Methods). We, next, focused on colorectal tumors to investigate which genes underlie the less connectivity trend. The genes in the paired normal networks of different colorectal tumor datasets, based on positive or negative coexpression, had mostly low overlap, shown by Jaccard index (Table 2). Only the positive and negative correlation based networks of the same dataset (such as Ade1Pos and Ade1Neg) showed a high overlap. Therefore, the networks which all lost connectivity, mostly shared a low fraction of members. Next, we analyzed the genes which had less degree value in tumor specific network compared to paired normal specific network. Jaccard index values of the different colorectal tumor networks as well as randomization based p-values for Jaccard index values were calculated (see Methods). There were mostly low overlap values between various differentially less connected genes, however, the overlap values were mostly statistically significant (Table 3). For instance, the positive and negative correlation based networks of Adenoma 2 dataset (Ade2Pos and Ade2Neg) had a Jaccard value of 0.85 with a p-value of 0.001, therefore the less connectivity trend of these two sample sets were mostly based on the same set of genes. However, the negative correlation based networks of Carcinoma 1 dataset (Car1Neg) and Adenoma 1 dataset (Ade1Neg) had a Jaccard value of 0.19 with a p-value of 0.965, therefore the less connectivity trend of these two sample sets were not based on a common set of genes. Hence, the presence of low Jaccard index values and high p-values showed that there was not a common set of gene list underlying the less connectivity trend in all colorectal tumor based networks. However, there were some common gene lists which could be underlying the less connectivity trend of some networks such as the various networks of Adenoma 1 and Adenoma 2 datasets.

Table 2

Jaccard index values of the coexpression networks of normal samples

Table 3

Jaccard index values of differentially less connected genes (upper value) and the p-values for the overlap (lower value)

The topology of the colorectal tumor gene coexpression networks were mostly weakly similar to scale-free networks (Table 1, see Appendix). Hub genes, the most highly connected genes, could be considered responsible for most of the connectivity of the networks. We examined the differential hub genes, which were the most highly differentially connected genes, thus greatly contributed to the loss of connectivity. Therefore, the less connectivity trend in tumor networks were investigated by analyzing the lists of differential hub genes. Top 5% of genes were selected as differential hubs, based on the difference of the degree values between tumor and normal networks. Differential hub gene lists were defined for two separate colorectal tumor datasets for positive and negative correlations. Pairwise comparison of differential hubs between different networks showed a very low overlap (Table 4). Jaccard index values ranged from 0.00 to a maximum of 0.13. However, some differential hub pairs had a significant overlap suggesting that some common hub genes could underlie the differential connectivity trend for some cases. However, there were also some nonsignificant hub pairs (Table 4). In addition, there was no common gene among the hub lists. There was not also any common gene among hubs, when positive and negative correlation based networks were considered separately. These values indicate that no common list of important differential hub genes could be defined, therefore, different genes were responsible for the less connectivity of different tumor networks. Furthermore, we selected the direct neighbors of differential hub genes to generate the hub neighborhoods. Pairwise comparison of the hub neighborhoods showed a greater overlap, with Jaccard index values ranging from 0.04 to 0.99 (Table 5). However, the overlap of hub neighborhoods was mostly nonsignificant, so that these overlap values were randomly expected (Table 5). Therefore, differential hubs mostly have different differential connections among various colorectal tumor based networks. Thus, the presence of low overlap values and the lack of all significant Jaccard index values indicate that no common gene neighborhood is responsible for the loss of connectivity in colorectal tumors.

Table 4

Jaccard index values of the differential hubs (top 5% most highly differentially less connected genes) (upper value) and the p-values for the overlap (lower value)

Table 5

Jaccard index values of the differential hub neighborhoods in the normal specific networks and the p-values for the overlap (lower value)

We also performed a functional gene set enrichment analysis for differential hub genes (see Methods). Interestingly, the differential hubs obtained from positive correlation based networks did not have any significant gene ontology terms, whereas the negative correlation based differential hubs had several terms such as ‘extracellular space’, ‘extracellular matrix’, etc. (Table 6). Ade1Neg, Ade2Neg and Car2Neg shared ‘extracellular exosome’ as a significant term. In addition, Car1Neg also had ‘extracellular exosome’ with 21% coverage, however its Benjamini corrected p-value was 0.08, therefore it was omitted in Table 6. Thus, all negative correlation based differential hubs had an association with extracellular exosome. When we analyzed the overlaps among the extracellular exosome associated negative correlation based differential hubs, there was no common gene among all 4 datasets; adenoma datasets did not share any gene, whereas carcinoma samples shared only 5 genes; FBXO2 (F-box protein 2), GPC1 (glypican 1), NEU1 (neuraminidase 1), PHGDH (phosphoglycerate dehydrogenase), and RARRES2 (retinoic acid receptor responder 2). Overall, differential hubs generated from negative correlation based networks had an enrichment of genes which were related to extracellular exosome, however, no common gene list was responsible for thisenrichment.

Table 6

Functional Gene Set Enrichment Values of Differential Hubs in Colorectal Samples and Common Methylation Targets of Carcinoma Networks

	Gene Set	Percentage of Genes (%)	Benjamini Corrected p-value
Ade 1 Pos	–	–	–
Ade 2 Pos	–	–	–
Ade 1 Neg	O-glycan processing	7.4	0.0003
	Extracellular space	22.3	0.0003
	Extracellular exosome	27.7	0.046
Ade 2 Neg	Extracellular matrix	8.0	0.0003
	Extracellular space	15.9	0.003
	Extracellular exosome	24.4	0.01
	Plasma membrane	30.1	0.017
	Heparin binding	6.2	0.0004
Car 1 Pos	–	–	–
Car 2 Pos	–	–	–
Car 1 Neg	Trans-Golgi network	3.4	0.025
Car 2 Neg	Oxidation-reduction process	10.7	0.0013
	Xenobiotic metabolic process	3.9	0.0075
	Platelet degranulation	3.9	0.048
	Extracellular exosome	26.2	0.0004
	Blood microparticle	4.4	0.0025
	Extracellular space	14.6	0.0056
	Pyridoxal phosphate binding	3.4	0.01
Common Methylation Targets of Carcinoma Networks	Extracellular matrix organization	9.2	0.001
	Planar cell polarity pathway involved in neural tube closure	3.3	0.037
	Digestive tract morphogenesis	3.3	0.032
	Extracellular space	20.8	0.0007
	Extracellular matrix	9.2	0.0019
	Basement membrane	5.0	0.0091
	Plasma membrane	36.7	0.016
	Sequence-specific DNA binding	14.2	0.0001
	Calcium ion binding	12.5	0.035
	RNA polymerase II core promoter proximal region sequence-specific DNA binding	8.3	0.049
	Heparin binding	5.8	0.043

Mutation and methylation targets with less connectivity in tumor

We analyzed the distribution of the degree values of mutation and methylation target genes for colorectal tumors. Mutation targets did not show a statistically significant degree distribution (see Appendix). Next, the degree distribution of methylation targets was compared to random gene lists. Methylation was previously linked to loss of connectivity in colorectal cancer [9]. Genes with significantly higher methylation levels in colon cancer were obtained from COSMIC database [10]. Methylation targets, showed a decrease of connectivity in tumor samples, similar to the overall network (Figs. 4, 5). When compared to random distribution, methylation targets, had higher connectivity in tumor samples of positive correlation based networks of colorectal carcinoma datasets (Figs. 4, 5). The connection and unique neighbor number distributions were close to being significant in the tumor networks of carcinoma 1, and carcinoma 2. There were 199 methylation targets in the colorectal carcinoma 1 tumor network, whereas there were 233 methylation targets in the colorectal carcinoma 2 tumor network. 120 of them were common. Functional gene set enrichment of these common methylation targets showed extracellular-related terms (Table 6). This was similar to the enrichment results for the negative correlation based differential hubs. There was no consistent significant distribution in the other networks (normal and tumor networks of adenoma 1 and 2, and normal networks of carcinoma 1 and 2). There was also no consistent significant distribution of methylation targets in negative correlation based networks (results not shown).

Fig.4

Total number of connections and unique number of neighbors of methylation targets compared to random lists for Normal (N) and Tumor (T) samples of positive correlation based networks of colorectal adenoma 1 (Ade 1; GSE8671 dataset) and colorectal carcinoma 1 (Car 1; GSE18105 dataset). The percentage of random values with lower values than methylation targets was shown.

Fig.5

Total number of connections and unique number of neighbors of methylation targets compared to random lists for Normal (N) and Tumor (T) samples of positive correlation based networks of colorectal adenoma 2 (Ade 2; GSE89076 dataset) and colorectal carcinoma 2 (Car 2; GSE89076 dataset). The percentage of random values with lower values than methylation targets was shown.

Conclusion

In this study, we showed a trend for lower connectivity between genes in the gene coexpression networks of colorectal adenoma and carcinoma as well as aldosterone producing adenoma and uterine leiomyoma. The results suggest a loss of coexpression in tumors, starting from precancerous stages. Our results confirm other studies suggesting that biological systems lose coherence when exposed to perturbations [11]. Thus, the loss of coherence starts from earlier stages than malignancy. We limited our analysis only to those datasets which included paired normal-tumor samples from patients, therefore our conclusions about the transformation of normal networks to tumor networks are stronger than conclusions made from nonpaired normal-tumor network analysis. We analyzed 4 different datasets of 3 different types of pre-cancerous stages, which was limited by available public data. Future work should include more pre-malignant tumors so that a stronger conclusion could be made for the system level connectivity change of different tumor types. Additionally, future studies could include more tumor stages to fully describe the stage by stage transformation of tumor networks, so that we could have a topological view of carcinogenesis. We also showed that no particular genes could be defined as being responsible for the loss of connectivity trend. There was mostly low or nonsignificant overlap among differentially less connected genes or differential hubs. Functional gene set enrichment analysis suggested exosome related genes to be enriched among differential hubs which were highly less connected in negative correlation based networks of tumor, when compared to paired normal networks. Lastly, methylation targets had higher connectivity values in the tumor networks of colorectal carcinoma samples. These genes also had enrichment of extracellular-associated gene ontology terms. Future studies could investigate more on the reasons for the loss of negative correlations of exosome related genes during tumorigenesis and the unexpectedly higher connectivity levels of extracellular space related methylation targets in colorectal carcinoma. Results from these studies could enable a better system level understanding of cancer so that more efficient perturbation strategies could be developed based on cancer system design principles.

Methods

Gene expression datasets

We collected patient-derived, paired tumor and normal (adjacent tissue to tumor) sample gene expression data from NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/). The datasets in this study included 2 datasets for colorectal adenoma (GSE8671, GSE89076) and 2 datasets for colorectal carcinoma (GSE18105, GSE89076), 1 dataset for aldosterone producing adenoma (GSE60042) and 1 dataset for uterine leiomyoma (GSE31699). For colorectal adenoma 1, GSM215051-GSM215114 of GSE8671 were used. For colorectal adenoma 2, GSM2358439, GSM2358440, GSM2358441, GSM2358442, GSM2358443, GSM2358444, GSM2358471, GSM2358472, GSM2358501, and GSM2358502 of GSE89076 were used. For colorectal carcinoma 1, GSM452629-GSM452662 of GSE18105 were used. For colorectal carcinoma 2, GSM2358438, GSM2358445-GSM2358470, GSM2358473-GSM2358516, GSM2358B7:BL7437 of GSE89076 were used. For aldosterone producing adenoma, GSM1464109- GSM1464122 of GSE60042 were used. For uterine leiomyoma, GSM786770- GSM786801 of GSE31699 were used. The datasets contained paired samples from the same patient such that direct comparison can be made between the tumor and normal samples. The platforms of GSE8671 and GSE18105 were Affymetrix Human Genome U133 Plus 2.0 Array. The platform of GSE89076 and GSE60042 were Agilent SurePrint G3 Human GE 8x60K Microarray. The platform of GSE31699 was Illumina HumanHT-12 V3.0 expression beadchip. The datasets contained paired samples from the same patient such that direct comparison can be made between the tumor and normal samples. For the Affymetrix Microarray based datasets (GSE8671, GSE18105), the raw values were obtained and normalized with Bioconductor Affy package Robust Multichip Average (RMA) method [12, 13]. For GSE89076 adenoma and carcinoma datasets, background subtracted processed signal intensity values were obtained and were converted to log values and quantile normalized. GSE8671 adenoma was denoted as adenoma 1, GSE18105 carcinoma was denoted as carcinoma 1, GSE89076 adenoma was denoted as adenoma 2 and GSE89076 carcinoma was denoted as carcinoma 2. For GSE31699 and GSE40062, already log transformed and quantile normalized values were obtained from NCBI GEO database.

There were 3 criteria for the pre-processing of the normalized datasets, before calculating pairwise correlations between genes. Firstly, there were multiple probeset values for a single gene (based on the NCBI Entrez Gene ID) in microarray datasets, from which we selected only one probeset with the maximum variation. Secondly, we removed the genes with average expression levels lower than a certain threshold (5), as very low values would cause false-positive correlations. Thirdly, we removed the genes with standard deviation values lower than a certain threshold (0.5), as very lowly varied values would also cause false-positive correlations. The arbitrary thresholds have no effect on our results as shown in the Appendix.

Co-expression networks

To construct networks with the most informative lists of genes, we collected highly and variably expressed genes (see Appendix). Using these genes, we constructed specific networks for paired normal and tumor samples from adenoma, leiomyoma and carcinoma datasets. To assess for co-expression, PCC values were calculated for normal and tumor samples separately. To test for the significance of the correlation, normal and tumor samples of a dataset were randomized and PCC values were re-calculated and a random value distribution was obtained. From this distribution a one-sided p-value was calculated for higher (or equal) values than the real (non-randomized) dataset PCC value. Then, the p-values were corrected for multiple testing by False Discovery Rate (FDR). FDR corrected p-value threshold of 0.05 was used to define significantly correlated gene pairs in normal or tumor samples separately. To denote a connection between a significantly correlated gene pair in tumor samples a differential correlation threshold was used; PCC value above a threshold (0.5) in tumor samples and a PCC value below a threshold (0.2) in paired normal samples were required. The same approach was taken to denote a connection in normal samples. Defining the connections in this way provided us with differential connections. Changing the preprocessing threshold values (for defining highly and variably expressed genes), the correlation threshold values, and even removing the differential correlation requirement had no effect on our results (see Appendix). Positive and negative gene-gene correlations were considered separately. We constructed weighted positive and negative; normal and tumor networks. In each network, we calculated the degree value of a gene, as the sum of significant correlations of a gene with other genes, excluding the self-correlation. Visualization and topological analysis of the networks were performed by using Cytoscape 3.7.1 [14].

Differentially less connected genes, differential hubs and mutation/methylation targets

Differential degree values of genes were defined based on the difference of the normal network degree value from the tumor network degree value. The genes with positive differential degree values were defined as differentially less connected genes. For selecting differential hubs, differentially less connected genes were ranked from the one with the highest degree difference to the one with lowest, top 5% genes from the list were defined as differential hubs. Overlap of network members, differentially less connected genes and differential hubs were analyzed by Jaccard index, which is the ratio of intersection to union for two sets. For significance testing, a background was selected (all genes in the normal network for the analysis of differentially less connected genes, all differentially less connected genes for the differential hubs) and random lists of the same size were generated and Jaccard index values for random sets were calculated, thus, a random value distribution was obtained. From this distribution a one-sided p-value was calculated for higher (or equal) values than the real (non-randomized) Jaccard index value. Lastly, p-values were corrected for multiple testing by FDR.

Functional gene set enrichment analysis of differential hubs was done by using DAVID 6.8 (https://david.ncifcrf.gov/home.jsp) [15, 16]. For colorectal adenoma 1 and carcinoma 1 datasets, Affymetrix Human Genome U133 Plus 2.0 Array, for colorectal adenoma 2 and carcinoma 2 datasets, Agilent HumanGenome background were used. Significant Gene Ontology terms (Benjamini p-value <  0.05) were shown in Table 6.

Mutation targets were obtained from Cancer Gene Census as of April 2018 (https://cancer.sanger.ac.uk/census) [17]. Both somatic and germline mutation targets of colorectal cancer were selected, which was a total of 50 genes. 11–25 out of 50 genes were in the differential correlation networks. Methylation targets were obtained from COSMIC database (https://cancer.sanger.ac.uk/cosmic) [10]. Genes with significantly higher levels of methylation in colon cancer, compared to normal, were obtained. We omitted the genes, which were grouped in both lower and higher methylation lists. Total number of methylation targets were 746. In the differential correlation networks 69–233 of them were present. The degree distribution of mutation targets was compared to random gene lists of the same size.

Author contributions

ED, OK, ZSO and CC designed the study, ED performed the computational analysis.