Building an infrastructure to support the development, conduct, and reporting of informative clinical studies: The Rockefeller University experience.

Introduction: Clinical trials are a vital component of translational science, providing crucial information on the efficacy and safety of new interventions and forming the basis for regulatory approval and/or clinical adoption. At the same time, they are complex to design, conduct, monitor, and report successfully. Concerns over the last two decades about the quality of the design and the lack of completion and reporting of clinical trials, characterized as a lack of "informativeness," highlighted by the experience during the COVID-19 pandemic, have led to several initiatives to address the serious shortcomings of the United States clinical research enterprise. Methods and Results: Against this background, we detail the policies, procedures, and programs that we have developed in The Rockefeller University Center for Clinical and Translational Science (CCTS), supported by a Clinical and Translational Science Award (CTSA) program grant since 2006, to support the development, conduct, and reporting of informative clinical studies. Conclusions: We have focused on building a data-driven infrastructure to both assist individual investigators and bring translational science to each element of the clinical investigation process, with the goal of both generating new knowledge and accelerating the uptake of that knowledge into practice.

Prediction of Liver Transplant Rejection With a Biologically Relevant Gene Expression Signature.

BACKGROUND: Noninvasive biomarkers distinguishing early immune activation before acute rejection (AR) could more objectively inform immunosuppression management in liver transplant recipients (LTRs). We previously reported a genomic profile distinguishing LTR with AR versus stable graft function. This current study includes key phenotypes with other causes of graft dysfunction and uses a novel random forest approach to augment the specificity of predicting and diagnosing AR. METHODS: Gene expression results in LTRs with AR versus non-AR (combination of other causes of graft dysfunction and normal function) were analyzed from single and multicenter cohorts. A 70:30 approach (61 ARs; 162 non-ARs) was used for training and testing sets. Microarray data were normalized using a LT-specific vector. RESULTS: Random forest modeling on the training set generated a 59-probe classifier distinguishing AR versus non-AR (area under the curve 0.83; accuracy 0.78, sensitivity 0.70, specificity 0.81, positive predictive value 0.54, negative predictive value [NPV] 0.89; F-score 0.61). Using a locked threshold, the classifier performed well on the testing set (accuracy 0.72, sensitivity 0.67, specificity 0.73, positive predictive value 0.48, NPV 0.86; F-score 0.56). Probability scores increased in samples preceding AR versus non-AR, when liver function tests were normal, and decreased following AR treatment (P < 0.001). Ingenuity pathway analysis of the genes revealed a high percentage related to immune responses and liver injury. CONCLUSIONS: We have developed a blood-based biologically relevant biomarker that can be detected before AR-associated graft injury distinct from LTR never developing AR. Given its high NPV ("rule out AR"), the biomarker has the potential to inform precision-guided immunosuppression minimization in LTRs.

Donor-derived cell-free DNA levels predict graft injury in liver transplant recipients.

Donor-derived cell-free DNA (dd-cfDNA) has been evaluated as a rejection marker in organ transplantation. This study sought to assess the utility of dd-cfDNA to diagnose graft injury in liver transplant recipients (LTR) and as a predictive biomarker prior to different causes of graft dysfunction. Plasma from single and multicenter LTR cohorts was analyzed for dd-cfDNA. Phenotypes of treated biopsy-proven acute rejection (AR, N = 57), normal function (TX, N = 94), and acute dysfunction no rejection (ADNR; N = 68) were divided into training and test sets. In the training set, dd-cfDNA was significantly different between AR versus TX (AUC 0.95, 5.3% cutoff) and AR versus ADNR (AUC 0.71, 20.4% cutoff). Using these cutoffs in the test set, the accuracy and NPV were 87% and 100% (AR vs. TX) and 66.7% and 87.8% (AR vs. ADNR). Blood samples collected serially from LTR demonstrated incremental elevations in dd-cfDNA prior to the onset of graft dysfunction (AR > ADNR), but not in TX. Dd-cfDNA also decreased following treatment of rejection. In conclusion, the serial elevation of dd-cfDNA identifies pre-clinical graft injury in the context of normal liver function tests and is greatest in rejection. This biomarker may help detect early signs of graft injury and rejection to inform LTR management strategies.

Epigenetic reprogramming of host and viral genes by Human Cytomegalovirus infection in Kasumi-3 myeloid progenitor cells at early times post-infection.

HCMV establishes latency in myeloid cells. Using the Kasumi-3 latency model, we previously showed that lytic gene expression is activated prior to establishment of latency in these cells. The early events in infection may have a critical role in shaping establishment of latency. Here, we have used an integrative multi-omics approach to investigate dynamic changes in host and HCMV gene expression and epigenomes at early times post infection. Our results show dynamic changes in viral gene expression and viral chromatin. Analyses of Pol II, H3K27Ac and H3K27me3 occupancy of the viral genome showed that 1) Pol II occupancy was highest at the MIEP at 4 hours post infection. However, it was observed throughout the genome; 2) At 24 hours, H3K27Ac was localized to the major immediate early promoter/enhancer and to a possible second enhancer in the origin of replication OriLyt; 3) viral chromatin was broadly accessible at 24 hpi. In addition, although HCMV infection activated expression of some host genes, we observed an overall loss of de novo transcription. This was associated with loss of promoter-proximal Pol II and H3K27Ac, but not with changes in chromatin accessibility or a switch in modification of H3K27.Importance.HCMV is an important human pathogen in immunocompromised hosts and developing fetuses. Current anti-viral therapies are limited by toxicity and emergence of resistant strains. Our studies highlight emerging concepts that challenge current paradigms of regulation of HCMV gene expression in myeloid cells. In addition, our studies show that HCMV has a profound effect on de novo transcription and the cellular epigenome. These results may have implications for mechanisms of viral pathogenesis.

A first-in-human phase 0 clinical study of RNA interference-based spherical nucleic acids in patients with recurrent glioblastoma.

Glioblastoma (GBM) is one of the most difficult cancers to effectively treat, in part because of the lack of precision therapies and limited therapeutic access to intracranial tumor sites due to the presence of the blood-brain and blood-tumor barriers. We have developed a precision medicine approach for GBM treatment that involves the use of brain-penetrant RNA interference-based spherical nucleic acids (SNAs), which consist of gold nanoparticle cores covalently conjugated with radially oriented and densely packed small interfering RNA (siRNA) oligonucleotides. On the basis of previous preclinical evaluation, we conducted toxicology and toxicokinetic studies in nonhuman primates and a single-arm, open-label phase 0 first-in-human trial (NCT03020017) to determine safety, pharmacokinetics, intratumoral accumulation and gene-suppressive activity of systemically administered SNAs carrying siRNA specific for the GBM oncogene Bcl2Like12 (Bcl2L12). Patients with recurrent GBM were treated with intravenous administration of siBcl2L12-SNAs (drug moniker: NU-0129), at a dose corresponding to 1/50th of the no-observed-adverse-event level, followed by tumor resection. Safety assessment revealed no grade 4 or 5 treatment-related toxicities. Inductively coupled plasma mass spectrometry, x-ray fluorescence microscopy, and silver staining of resected GBM tissue demonstrated that intravenously administered SNAs reached patient tumors, with gold enrichment observed in the tumor-associated endothelium, macrophages, and tumor cells. NU-0129 uptake into glioma cells correlated with a reduction in tumor-associated Bcl2L12 protein expression, as indicated by comparison of matched primary tumor and NU-0129-treated recurrent tumor. Our results establish SNA nanoconjugates as a potential brain-penetrant precision medicine approach for the systemic treatment of GBM.

Identification of Geographic Specific SARS-Cov-2 Mutations by Random Forest Classification and Variable Selection Methods.

RNA viral genomes have very high mutations rates. As infection spreads in the host populations, different viral lineages emerge acquiring independent mutations that can lead to varied infection and death rates in different parts of the world. By application of Random Forest classification and feature selection methods, we developed an analysis pipeline for identification of geographic specific mutations and classification of different viral lineages, focusing on the missense-variants that alter the function of the encoded proteins. We applied the pipeline on publicly available SARS-CoV-2 datasets and demonstrated that the analysis pipeline accurately identified country or region-specific viral lineages and specific mutations that discriminate different lineages. The results presented here can help designing country-specific diagnostic strategies and prioritizing the mutations for functional interpretation and experimental validations.

FTO-Dependent N 6-Methyladenosine Modifications Inhibit Ovarian Cancer Stem Cell Self-Renewal by Blocking cAMP Signaling.

N 6-Methyladenosine (m6A) is the most abundant modification of mammalian mRNAs. RNA methylation fine tunes RNA stability and translation, altering cell fate. The fat mass- and obesity-associated protein (FTO) is an m6A demethylase with oncogenic properties in leukemia. Here, we show that FTO expression is suppressed in ovarian tumors and cancer stem cells (CSC). FTO inhibited the self-renewal of ovarian CSC and suppressed tumorigenesis in vivo, both of which required FTO demethylase activity. Integrative RNA sequencing and m6A mapping analysis revealed significant transcriptomic changes associated with FTO overexpression and m6A loss involving stem cell signaling, RNA transcription, and mRNA splicing pathways. By reducing m6A levels at the 3'UTR and the mRNA stability of two phosphodiesterase genes (PDE1C and PDE4B), FTO augmented second messenger 3', 5'-cyclic adenosine monophosphate (cAMP) signaling and suppressed stemness features of ovarian cancer cells. Our results reveal a previously unappreciated tumor suppressor function of FTO in ovarian CSC mediated through inhibition of cAMP signaling. SIGNIFICANCE: A new tumor suppressor function of the RNA demethylase FTO implicates m6A RNA modifications in the regulation of cyclic AMP signaling involved in stemness and tumor initiation.

Discovery and validation of a novel blood-based molecular biomarker of rejection following liver transplantation.

Noninvasive biomarker profiles of acute rejection (AR) could affect the management of liver transplant (LT) recipients. Peripheral blood was collected following LT for discovery (Northwestern University [NU]) and validation (National Institute of Allergy and Infectious Diseases Clinical Trials in Organ Transplantation [CTOT]-14 study). Blood gene profiling was paired with biopsies showing AR or ADNR (acute dysfunction no rejection) as well as stable graft function samples (Transplant eXcellent-TX). CTOT-14 subjects had serial collections prior to AR, ADNR, TX, and after AR treatment. NU cohort gene expression (46 AR, 45 TX) was analyzed using random forest models to generate a classifier training set (36 gene probe) distinguishing AR vs TX (area under the curve 0.92). The algorithm and threshold were locked and tested on the CTOT-14 validation cohort (14 AR, 50 TX), yielding an accuracy of 0.77, sensitivity 0.57, specificity 0.82, positive predictive value (PPV) 0.47, and negative predictive value (NPV) 0.87 for AR vs TX. The probability score line slopes were positive preceding AR, and negative preceding TX and non-AR (TX + ADNR) (P ≤ .001) and following AR treatment. In conclusion, we have developed a blood biomarker diagnostic for AR that can be detected prior to AR-associated graft injury as well a normal graft function (non-AR). Further studies are needed to evaluate its utility in precision-guided immunosuppression optimization following LT.

Immune landscapes associated with different glioblastoma molecular subtypes.

Recent work has highlighted the tumor microenvironment as a central player in cancer. In particular, interactions between tumor and immune cells may help drive the development of brain tumors such as glioblastoma multiforme (GBM). Despite significant research into the molecular classification of glioblastoma, few studies have characterized in a comprehensive manner the immune infiltrate in situ and within different GBM subtypes.In this study, we use an unbiased, automated immunohistochemistry-based approach to determine the immune phenotype of the four GBM subtypes (classical, mesenchymal, neural and proneural) in a cohort of 98 patients. Tissue Micro Arrays (TMA) were stained for CD20 (B lymphocytes), CD5, CD3, CD4, CD8 (T lymphocytes), CD68 (microglia), and CD163 (bone marrow derived macrophages) antibodies. Using automated image analysis, the percentage of each immune population was calculated with respect to the total tumor cells. Mesenchymal GBMs displayed the highest percentage of microglia, macrophage, and lymphocyte infiltration. CD68+ and CD163+ cells were the most abundant cell populations in all four GBM subtypes, and a higher percentage of CD163+ cells was associated with a worse prognosis. We also compared our results to the relative composition of immune cell type infiltration (using RNA-seq data) across TCGA GBM tumors and validated our results obtained with immunohistochemistry with an external cohort and a different method. The results of this study offer a comprehensive analysis of the distribution and the infiltration of the immune components across the four commonly described GBM subgroups, setting the basis for a more detailed patient classification and new insights that may be used to better apply or design immunotherapies for GBM.

Transcriptional consequences of impaired immune cell responses induced by cystic fibrosis plasma characterized via dual RNA sequencing.

BACKGROUND: In cystic fibrosis (CF), impaired immune cell responses, driven by the dysfunctional CF transmembrane conductance regulator (CFTR) gene, may determine the disease severity but clinical heterogeneity remains a major therapeutic challenge. The characterization of molecular mechanisms underlying impaired immune responses in CF may reveal novel targets with therapeutic potential. Therefore, we utilized simultaneous RNA sequencing targeted at identifying differentially expressed genes, transcripts, and miRNAs that characterize impaired immune responses triggered by CF and its phenotypes. METHODS: Peripheral blood mononuclear cells (PBMCs) extracted from a healthy donor were stimulated with plasma from CF patients (n = 9) and healthy controls (n = 3). The PBMCs were cultured (1 × 105 cells/well) for 9 h at 37 ° C in 5% CO2. After culture, total RNA was extracted from each sample and used for simultaneous total RNA and miRNA sequencing. RESULTS: Analysis of expression signatures from peripheral blood mononuclear cells induced by plasma of CF patients and healthy controls identified 151 genes, 154 individual transcripts, and 41 miRNAs differentially expressed in CF compared to HC while the expression signatures of 285 genes, 241 individual transcripts, and seven miRNAs differed due to CF phenotypes. Top immune pathways influenced by CF included agranulocyte adhesion, diapedesis signaling, and IL17 signaling, while those influenced by CF phenotypes included natural killer cell signaling and PI3K signaling in B lymphocytes. Upstream regulator analysis indicated dysregulation of CCL5, NF-κB and IL1A due to CF while dysregulation of TREM1 and TP53 regulators were associated with CF phenotype. Five miRNAs showed inverse expression patterns with three target genes relevant in CF-associated impaired immune pathways while two miRNAs showed inverse expression patterns with two target genes relevant to a dysregulated immune pathway associated with CF phenotypes. CONCLUSIONS: Our results indicate that miRNAs and individual transcript variants are relevant molecular targets contributing to impaired immune cell responses in CF.

Platform-Independent Classification System to Predict Molecular Subtypes of High-Grade Serous Ovarian Carcinoma.

PURPOSE: Molecular cancer subtyping is an important tool in predicting prognosis and developing novel precision medicine approaches. We developed a novel platform-independent gene expression-based classification system for molecular subtyping of patients with high-grade serous ovarian carcinoma (HGSOC). METHODS: Unprocessed exon array (569 tumor and nine normal) and RNA sequencing (RNA-seq; 376 tumor) HGSOC data sets, with clinical annotations, were downloaded from the Genomic Data Commons portal. Sample clustering was performed by non-negative matrix factorization by using isoform-level expression estimates. The association between the subtypes and overall survival was evaluated by Cox proportional hazards regression model after adjusting for the covariates. A novel classification system was developed for HGSOC molecular subtyping. Robustness and generalizability of the gene signatures were validated using independent microarray and RNA-seq data sets. RESULTS: Sample clustering recaptured the four known The Cancer Genome Atlas molecular subtypes but switched the subtype for 22% of the cases, which resulted in significant (P = .006) survival differences among the refined subgroups. After adjusting for covariate effects, the mesenchymal subgroup was found to be at an increased hazard for death compared with the immunoreactive subgroup. Both gene- and isoform-level signatures achieved more than 92% prediction accuracy when tested on independent samples profiled on the exon array platform. When the classifier was applied to RNA-seq data, the subtyping calls agreed with the predictions made from exon array data for 95% of the 279 samples profiled by both platforms. CONCLUSION: Isoform-level expression analysis successfully stratifies patients with HGSOC into groups with differing prognosis and has led to the development of robust, platform-independent gene signatures for HGSOC molecular subtyping. The association of the refined The Cancer Genome Atlas HGSOC subtypes with overall survival, independent of covariates, enhances the clinical annotation of the HGSOC cohort.

A clinically relevant murine model unmasks a "two-hit" mechanism for reactivation and dissemination of cytomegalovirus after kidney transplant.

Reactivation of latent cytomegalovirus remains an important complication after transplant. Although immunosuppression (IS) has been implicated as a primary cause, we have previously shown that the implantation response of a kidney allograft can lead to early transcriptional activation of latent murine cytomegalovirus (MCMV) genes in an immune-competent host and to MCMV reactivation and dissemination to other organs in a genetically immune-deficient recipient. We now describe a model that allows us to separately analyze the impact of the implantation effect vs that of a clinically relevant IS regimen. Treatment with IS of latently infected mice alone does not induce viral reactivation, but transplant of latently infected allogeneic kidneys combined with IS facilitates MCMV reactivation in the graft and dissemination to other organs. The IS regimen effectively dampens allo-immune inflammatory pathways and depletes recipient anti-MCMV but does not affect ischemia-reperfusion injury pathways. MCMV reactivation similar to that seen in allogeneic transplants combined with also occurs after syngeneic transplants. Thus, our data strongly suggest that while ischemia-reperfusion injury of the implanted graft is sufficient and necessary to initiate transcriptional reactivation of latent MCMV ("first hit"), IS is permissive to the first hit and facilitates dissemination to other organs ("second hit").

Development and clinical validity of a novel blood-based molecular biomarker for subclinical acute rejection following kidney transplant.

Noninvasive biomarkers are needed to monitor stable patients after kidney transplant (KT), because subclinical acute rejection (subAR), currently detectable only with surveillance biopsies, can lead to chronic rejection and graft loss. We conducted a multicenter study to develop a blood-based molecular biomarker for subAR using peripheral blood paired with surveillance biopsies and strict clinical phenotyping algorithms for discovery and validation. At a predefined threshold, 72% to 75% of KT recipients achieved a negative biomarker test correlating with the absence of subAR (negative predictive value: 78%-88%), while a positive test was obtained in 25% to 28% correlating with the presence of subAR (positive predictive value: 47%-61%). The clinical phenotype and biomarker independently and statistically correlated with a composite clinical endpoint (renal function, biopsy-proved acute rejection, ≥grade 2 interstitial fibrosis, and tubular atrophy), as well as with de novo donor-specific antibodies. We also found that <50% showed histologic improvement of subAR on follow-up biopsies despite treatment and that the biomarker could predict this outcome. Our data suggest that a blood-based biomarker that reduces the need for the indiscriminate use of invasive surveillance biopsies and that correlates with transplant outcomes could be used to monitor KT recipients with stable renal function, including after treatment for subAR, potentially improving KT outcomes.

Metastatic cancers promote cachexia through ZIP14 upregulation in skeletal muscle.

Patients with metastatic cancer experience a severe loss of skeletal muscle mass and function known as cachexia. Cachexia is associated with poor prognosis and accelerated death in patients with cancer, yet its underlying mechanisms remain poorly understood. Here, we identify the metal-ion transporter ZRT- and IRT-like protein 14 (ZIP14) as a critical mediator of cancer-induced cachexia. ZIP14 is upregulated in cachectic muscles of mice and in patients with metastatic cancer and can be induced by TNF-α and TGF-ß cytokines. Strikingly, germline ablation or muscle-specific depletion of Zip14 markedly reduces muscle atrophy in metastatic cancer models. We find that ZIP14-mediated zinc uptake in muscle progenitor cells represses the expression of MyoD and Mef2c and blocks muscle-cell differentiation. Importantly, ZIP14-mediated zinc accumulation in differentiated muscle cells induces myosin heavy chain loss. These results highlight a previously unrecognized role for altered zinc homeostasis in metastatic cancer-induced muscle wasting and implicate ZIP14 as a therapeutic target for its treatment.

Association of increased primary breast tumor AGR2 with decreased disease-specific survival.

OBJECTIVE: Tumor expression of Anterior Gradient 2 (AGR2), an endoplasmic reticulum protein disulfide isomerase, was associated with decreased breast cancer survival. We aimed to validate the association of tumor AGR2 mRNA expression with disease-specific survival (DSS) and identify differentially expressed signaling pathways between high and low AGR2 expression tumor groups. METHODS: Primary tumor mRNA expression data from the METABRIC study was used to evaluate AGR2 expression as a prognostic factor for DSS while adjusting for survival-determining confounders using Cox proportional-hazards regression. Differentially expressed genes and signaling pathway differences between high and low AGR2 groups were determined by modular enrichment analyses using DAVID and Ingenuity Pathway Analysis. RESULTS: Increased tumor AGR2 mRNA expression was associated with decreased DSS among 1,341 women (per each standard deviation increase of AGR2 expression: HR 1.14, 95% CI: 1.01-1.29, P = 0.03). Pathway analyses supported prior experimental studies showing that estrogen receptor 1 (ESR1) regulated AGR2 expression. Canonical signaling pathways significantly differentially represented between high and low AGR2 groups included those involved in inflammation and immunity. CONCLUSION: Increased primary tumor AGR2 expression was associated with decreased DSS. Pathway analyses suggested that increased AGR2 was associated with endoplasmic reticular homeostasis, possibly allowing tumor cells to overcome hypoxic stress and meet the increased protein demand of tumorigenesis, thereby preventing unfolded protein response-mediated apoptosis.

Comparative evaluation of isoform-level gene expression estimation algorithms for RNA-seq and exon-array platforms.

Given that the majority of multi-exon genes generate diverse functional products, it is important to evaluate expression at the isoform level. Previous studies have demonstrated strong gene-level correlations between RNA sequencing (RNA-seq) and microarray platforms, but have not studied their concordance at the isoform level. We performed transcript abundance estimation on raw RNA-seq and exon-array expression profiles available for common glioblastoma multiforme samples from The Cancer Genome Atlas using different analysis pipelines, and compared both the isoform- and gene-level expression estimates between programs and platforms. The results showed better concordance between RNA-seq/exon-array and reverse transcription-quantitative polymerase chain reaction (RT-qPCR) platforms for fold change estimates than for raw abundance estimates, suggesting that fold change normalization against a control is an important step for integrating expression data across platforms. Based on RT-qPCR validations, eXpress and Multi-Mapping Bayesian Gene eXpression (MMBGX) programs achieved the best performance for RNA-seq and exon-array platforms, respectively, for deriving the isoform-level fold change values. While eXpress achieved the highest correlation with the RT-qPCR and exon-array (MMBGX) results overall, RSEM was more highly correlated with MMBGX for the subset of transcripts that are highly variable across the samples. eXpress appears to be most successful in discriminating lowly expressed transcripts, but IsoformEx and RSEM correlate more strongly with MMBGX for highly expressed transcripts. The results also reinforce how potentially important isoform-level expression changes can be masked by gene-level estimates, and demonstrate that exon arrays yield comparable results to RNA-seq for evaluating isoform-level expression changes.

Fluorescence and piezoresistive cantilever sensing of trinitrotoluene by an upper-rim tetrabenzimidazole conjugate of calix[4]arene and delineation of the features of the complex by molecular dynamics.

A new benzimidazole-functionalized calix[4]arene receptor (R) was synthesized and characterized. The receptor R shows better selectivity toward trinitrotoluene (TNT) compared to the other nitro explosives in solution, which also retains its effectiveness for solid-phase detection. The chemical interactions of the molecule with different nitro explosive analytes were studied by fluorescence spectroscopy and by a molecular dynamics approach. The molecular dynamics studies show a 1:3 complex between R and TNT, and hence high sensitivity was imparted by fluorescence studies. The detection of explosive vapors in ambient conditions was tested by using a sensitive coating layer of R on an SU-8/CB-based piezoresistive cantilever surface. The developed device showed large sensitivity toward TNT compared to cyclotrimethylenetrinitramine (RDX) and pentaerythritol tetranitrate (PETN) in the solid state at their respective vapor pressures at room temperature. The detection sensitivity of the device was estimated to be 35 mV for TNT at ambient conditions. Moreover, the sensor does not show a response when exposed to humidity. These results demonstrate that R can be used as one of the coating materials for a cantilever for the detection of TNT using piezoresistivity measurement. R can also detect the explosives in solution with high sensitivity and selectivity by fluorescence spectroscopy.

Genetic programming-based approach to elucidate biochemical interaction networks from data.

Biochemical systems are characterised by cyclic/reversible reciprocal actions, non-linear interactions and a mixed relationship structures (linear and non-linear; static and dynamic). Deciphering the architecture of such systems using measured data to provide quantitative information regarding the nature of relationships that exist between the measured variables is a challenging proposition. Causality detection is one of the methodologies that are applied to elucidate biochemical networks from such data. Autoregressive-based modelling approach such as granger causality, partial directed coherence, directed transfer function and canonical variate analysis have been applied on different systems for deciphering such interactions, but with limited success. In this study, the authors propose a genetic programming-based causality detection (GPCD) methodology which blends evolutionary computation-based procedures along with parameter estimation methods to derive a mathematical model of the system. Application of the GPCD methodology on five data sets that contained the different challenges mentioned above indicated that GPCD performs better than the other methods in uncovering the exact structure with less false positives. On a glycolysis data set, GPCD was able to fill the 'interaction gaps' which were missed by other methods.