The task of learning a causal model from observational data, or a combination of observational and interventional data, is commonly referred to as a causal discovery or causal structure learning [1] . Causal discovery from two variables based on observational data in the absence of time series or controlled interventions is a challenging problem and necessitates additional assumptions [2] . This is a ubiquitous problem in almost all domains of science, but particularly so in econometrics, meteorology, biology and medicine where interventional approaches are difficult or in several cases not feasible.
Model-free data-driven approaches for causal discovery have developed significantly over the past decade or so in an attempt to address the problem of causal discovery such as Granger Causality (GC) [3] , Transfer Entropy (TE) [4] and Compression-Complexity Causality (CCC) [5] . These methods have been used in various disciplines across neuroscience, climatology, econometrics, etc and rely on properties of time-series data. Both GC and TE have assumptions that need to be met for satisfactory inference, while CCC is assumption-free and robust to many artefacts and nuisance variables. All three need careful parameter calibration and selection for optimally accurate performance.
A class of model-free causal discovery methods do not assume a temporal structure in the data and are rooted in algorithmic information theory, chiefly based on the notion of Kolmogorov complexity. The Kolmogorov complexity of a finite binary string is the length of the shortest binary program that generates that string and reflects the computational resources needed to specify it. For two observed variables X and Y , causal inference can be made by identifying the direction between X and Y where factorization of their joint distribution yields the lowest total Kolmogorov complexity [6] . Since Kolmogorov complexity is not computable, it has typically been approximated using the Minimum Description Length (MDL) principle [7] . MDL relies on lossless compressors to estimate the shortest program or description.
ORIGO is a method based on the MDL principle which relies on causal inference by tree-based compression [6] . ORIGO infers that X is a likely cause of Y if better compression is achieved by first compressing X and then compressing Y given X, than in the other direction. The crux of this approach relies on the postulate that if X causes Y , then describing Y using X will be easier or simpler than in the other direction [8] . In other words, if X causes Y , X will provide more information about Y than vice versa. This approach is employed in ERGO [9] which uses complexity estimates in lieu of Kolmogorov complexity instead of lossless compression. Both ORIGO and ERGO have been validated and benchmarked for univariate causal discovery.
Information theoretic quantities, such as Shannon entropy, suffer from several limitations for characterizing dynamical complexity of short and noisy time series data and better alternatives that capture complexity have been proposed and rigorously tested [10] . These are called Compression-Complexity Measures (CCMs) and make use of lossless compression algorithms to characterize complexity, and has been validated using Lempel-Ziv (LZ) complexity [11] and Effort-to-Compress (ETC) algorithms [12] . Both LZ and ETC are grammar-based compression algorithms which infer a context-free grammar (CFG) based on an input sequence [13] and have found diverse applications spanning numerous domains. These CCMs have been demonstrated to be robust to noise, artefacts, missing data and further shown to capture complex behavior of dynamical and stochastic systems [10] .
In the present paper, we propose a framework for causal discovery based on these CCMs that rely on inferred grammars to quantify information exchanged and shared between a pair of sequences. Based on this framework, we propose three models and evaluate their performance on synthetic as well as real-world benchmark data and compare their performance with that of ORIGO and ERGO.
Further, since our usage of CCMs within the proposed framework permits valid causal inference from sequences with different lengths, we propose and carry out a unique application of our models for assessing directional information exchange between genomic sequences of 16,619 SARS-CoV-2 virus isolates obtained from human samples. We believe this novel application has the potential to facilitate further investigations concerning key issues in bioinformatics, systems biology and epidemiology involving contact-tracing, epidemic monitoring, evolution and genetic interactions.
We first propose our framework for univariate causal discovery from discrete symbolic sequences, which consist of ordered sets of elements or symbols, recorded with or without a concrete notion of time, for example: {2, 3, 1, 2, 3, 2...}, {a, b, b, a, c, a, b, ...}. We then present models based on this framework followed by empirical applications to synthetic and real-world data. We also present a novel experimental application of this framework for genome sequences.
Grammar-based compression schemes construct a context-free grammar (CFG) for a given string x to be compressed by transforming x into a CFG, G [13] . The compressed representation of x itself has been used for extracting information from and searching in sequences. We focus on G, which is a generative model of x and contains its parsed structure based on rules relating to elements of x.
For two non-identical discrete symbolic sequences x and y with the same lengths and set of symbols, a lossless grammar-based compressor, L, will construct their CFGs, G x and G y respectively. G x represents an ideal generative grammar for x, and G y for y. L will output a compressed representation L(x|G x ) or L(x) of x and analogously L(y|G y ) or L(y) of y. Since x and y are distinct sequences, G x and G y are also distinct. Using G x to parse or compress y, L(y|G x ) leads to non-ideal compression of y compared to L(y). Consequently, the compressed representation of y can potentially be larger if compressed using G x instead of G y . Such an approach in the context of the Minimum Description Length (MDL) principle [7] has been employed for ORIGO [6] .
The compressed representations of x and y can be described alternatively in terms of Compression-Complexity Measures (CCMs) [10] besides description lengths. CCMs are measures of complexity derived from lossless data compression algorithms and have been shown to be robust estimators of the complexity of discrete as well as continuous dynamical ans stochastic systems. In the present context, we argue that the CCM of L(y|G x ) is different from the CCM of L(y). G x may compress y better or poorer than G y , resulting in CCM of L(y|G x ) lesser or greater than that of L(y) respectively. Similarly estimating CCMs of L(x|G y ) and L(x) and comparing both directions permits an examination of the influence of G x and G y on the sequences x and y. We extend this notion further and present 2 formulations for causal discovery, one based on penalty and another based on efficacy of compression using non-ideal grammars.
The difference CCM L (y|G x )−CCM L (y) represents the cost in terms of compression-complexity incurred by compressing y using a non-ideal grammar, in this case G x . The better G x can compress y, the smaller the CCM of L(y|G x ) and thus, smaller the cost. This cost or penalty, P , is directional and consists of the penalty of compressing y using G x , P x→y as well as the penalty of compressing x using G y , P y→x .
If penalty in one direction is lesser than that in the other direction, then one sequence's grammar can better compress the other sequence. This implies that the if P x→y < P y→x then the inferred generative model G x can account for y better than G y can account for x. We extend this to arrive at the following rules for inferring a causal direction:
If P x→y < P y→x , we infer x → y, If P x→y > P y→x , we infer y → x, If P x→y = P y→x , we are undecided.
A threshold can be introduced for the differences for inferring direction in practice. The difference |P y→x − P x→y | is an indication of the causal strength or strength of causal evidence in favor of the inferred direction, the larger the difference, the stronger the causal evidence in the inferred direction.
The difference CCM L (y|G x )−CCM L (x) represents the efficacy of G x towards complete lossless compression of y. If G x compresses y effectively, then residual y will have lower CCM and this difference will be small. The difference can also be viewed as the additional compression-complexity introduced by using G x to compress y instead of x. This efficacy, E is directional and consists of the efficacy of compressing y using G x , E x→y as well as the efficacy of compressing x using G y , E y→x .
If efficacy in one direction is greater than that in the other direction, then one sequence's grammar can more effectively compress the other sequence. This implies that the if E x→y > E y→x then the inferred generative model G x can account for y better than G y can account for x. We extend this to arrive at the following rules for inferring a causal direction:
As with the penalty approach, a threshold can be introduced for the differences for inferring direction here as well. The difference |E y→x − E x→y | indicates causal strength or strength of causal evidence in favor of the inferred direction, the larger the difference, the stronger the causal evidence in the inferred direction.
The framework presented involves describing a sequence and its compressed representation using CCMs. While various lossless grammar-based compression schemes exist, CCMs have been studied and validated using Lempel-Ziv (LZ) complexity and Effort-To-Compress (ETC) measures [10] . We present implementations of our framework using these two CCMs for both the penalty and efficacy approaches. We also use and describe the ORIGO and ERGO models from two similar frameworks rooted in algorithmic information theory, for comparative assessment of performance.
ETC is defined as the effort to compress [12] an input symbolic sequence using the lossless compression algorithm known as Non-Sequential Recursive Pair Substitution (NSRPS) [14] . Numerically ETC is an estimate of the number of steps required by NSRPS to compress an input sequence to a constant sequence or a sequence with zero entropy.
ETC has been demonstrated to reliably capture complexity of short and noisy time series and is robust to missing data [15, 10] . The algorithm is functionally similar to the Re-Pair lossless compressor [16, 17] although ETC has found broad applications in diverse domains including signal processing, de-noising, cognitive studies, heart-rate variability, Schizophrenia research, etc [18, 19, 20, 21, 22, 23, 24] . ETC has also been employed as the compression-complexity measure underlying the Compression-Complexity Causality (CCC) framework for robust causal inference through data-based intervention [5] .
Since ETC is an off-line compressor, G x as well as G y can be stored separately and residual y and residual x can be observed individually, simplifying the formulation to:
The efficacy of G x for compressing y and vice-versa can be assessed by directly estimating CCMs of the residual sequences, leading to the formulation:
where λ is a normalization factor to allow for better comparability of ETC estimates of the residuals.
We use the standard normalization factor for ETC where λ = (length(y residual ) − 1) −1 , referred to as the ETC-E formulation.
Lempel-Ziv complexity (1976) is defined as the number of different sub-strings encountered as the binary sequence is viewed as a stream (from left to right) [11] . LZ complexity has been used extensively in a wide spectrum of domains from linguistics, phylogenetics and neural-spike train analysis to mechanical fault identification and siezure detection [25, 26, 27, 28, 29] . We demonstrate the usage of LZ complexity in the penalty model, where CCM L (y|G x ) can be directly obtained by concatenation of the two sequences as xy. Here, LZ compresses the concatenated sequence as a stream and infers a grammar from x first, until it encounters y which is compressed using the existing grammar as well as an updated grammar when new substrings are encountered in y. This permits the formulation:
ORIGO is an efficient method for causal inference [6] from binary data based on the minimum description length (MDL) principle [7] . It relies on the tree-based lossless compressor, PACK [30] , for searching a set of decision trees that encode the data most succinctly. ORIGO works with both univariate and multivariate data with an equal number of observations and has been extensively evaluated previously. We chose ORIGO for comparison since it works directly on the data and does not require assumptions about distributions, similar to our proposed framework.
ERGO is a robust method for causal inference based on direction on information and uses cumulative and Shannon entropy [9] . It determines the amount of information one set of data provides about another and vice versa, and infers causal direction based on the strongest direction of information. ERGO also works with both univariate and multivariate data and has been shown to be robust to noise. We chose ERGO for comparison due to the principles of causal inference it shares with our proposed framework based on direction of information. While we employ robust CCMs to estimate strengths under the penalty and efficacy models, ERGO uses conditional cumulative entropy to estimate relative amount of information shared based on Kolmogorov Complexity.
An open-source implementation of ETC in Python (version 3.8) was used, with some steps of NSRPS implemented in Cython [31] and NumPy [32] , available on GitHub: https://github.com/pranaysy/ETCPy/ Simulations and processing of data used in the paper were carried out using Python scripts. LZ complexity was computed using the implementation provided in the open-source package EntroPy for Python on GitHub: https://github.com/raphaelvallat/entropy
Estimates for ORIGO and ERGO were computed using the reference open-source implementation in Python and Java respectively, available publicly at http://eda.mmci.uni-saarland.de/prj/origo/
The autoregressive (AR) model is widely used in statistics, econometrics as well as signal processing ans is used to describe time-varying processes that linearly depend on their own past as well as the past of other processes [33] . The AR model is commonly used for tests of causal inference and is foundational to the Granger causal framework [3] . We used the AR model for unidirectional causal inference and simulated autoregressive processes of order one (AR(1)) as follows with X and Y as the dependent and independent processes respectively.
where a = b = 0.8, t = 1 to 1,000s, sampling period = 1s. We varied the coupling parameter, φ, from 0 to 0.95 in steps of 0.05 for a total of 20 values. Noise terms X , Y = vN , where noise intensity, v = 0.01 and follows the standard normal distribution. For each value of φ, we randomly sampled 1,000 trials of X-Y pairs, and discretized each sampled X-Y pair to binary sequences using an equi-width binning strategy. We then estimated causal direction and scores using all five models for each trial, for a total of 20,000 trials.
We evaluated performance as accuracy against decision rates across varying coupling strengths for each model. This was calculated as accuracy over top-k% X-Y pairs after sorting in descending order the absolute difference in estimated scores from each of the five models. For pairs where a direction could not be discovered, we flip a fair coin. We also present overall accuracy based on 1,000 trials for each coupling strength for all five models along with estimates of area under the Receiver Operating Characteristic curve (AUROC) and area under the Precision-Recall curve (AUPRC).
We evaluated the performance of all five models on benchmark cause-effect pairs with known ground truth, which consists of real-world causal variables [34] . The dataset consists of both multivariate and univariate cause-effect pairs, and here we consider only the 90 univariate ones. In order to facilitate comparability, we used the interaction-preserving discretization (IPD) strategy [35] for all the pairs as used previously [6] . Similar to the evaluation of performance for synthetic data, we investigated accuracy against decision rates as well as AUROC and AUPRC for all five models.
The highly contagious pathogen called Severe Acute Respiratory syndrome coronavirus 2 (SARS-CoV-2) was first identified in December 2019 [36] as the cause of a respiratory illness designated coronavirus disease 2019, or Covid-19, classified as a global pandemic by the World Health Organization (WHO) on 11th March 2020. Coronaviruses are generally well known to evolve environmentally and have high virulence conferring effective transmission and immune evasion strategies [37] . At the time of writing (October 15th, 2020) over 38 million cases and 1.09 million fatalities have been reported globally [38] , with neither a cure for Covid-19 nor an effective vaccine within sight [39] .
Effective countermeasures against the virus require the development of data and tools to understand and monitor its spread and immune responses to it [40] . Methods grounded in information theory and data compression have found numerous applications in genomics -global sequence analysis, phylogenetics, evolutionary modeling, sequence complexity, motif discovery and classification, analysis of secondary structures, etc [41, 42] . Recently, approaches for causal inference that rely on graphical models [43] as well as additive noise models (ANMs) [44] have received attention for causal genetic analysis and bivariate causal discovery from gene expression data [45] . Here, we present a unique application of our proposed grammar-based framework for inferring a direction of causal information exchange directly from genome sequences obtained from viral isolates of SARS-CoV-2.
We hypothesized that the SARS-CoV-2 consensus sequence [46] would capture information representative of the initial state of the viral genome prior to the global outbreak. Subsequent transmission of the virus would be accompanied with genome evolution over time, altering the nucleotide sequence's information content. Since this evolved sequence contains changes accumulated over the consensus sequence, the evolved sequence can be viewed as a derivative of the consensus sequence, admitting a direction. Alternatively then, the consensus sequence potentially 'causes' the evolved sequence, permitting us to examine this within the causal discovery framework and test if our models indeed identify this hypothesized direction. Similarly, we also examined if the first sequence reported in each country 'caused' subsequent sequences in that country.
To evaluate this hypothesis, we computed estimates from the three proposed models of our framework -ETC-E, ETC-P and LZ-P. Since pairs of sequences are never evaluated simultaneously at a symbol-level under these three models, novel applications to genomic data are possible where sequences with different lengths are the norm. Pairwise estimates were computed for sequences from each country such that one sequence of a pair was first fixed to the consensus sequence, and later fixed to the first sequence for that country, resulting in estimates of direction and strength. Based on these estimates, we evaluated the hypothesis involving causal directions using the proportion of sequences for each country for which causal discovery was made in the expected direction for both the consensus sequence as well as the first sequence per country.
These experiments involved estimation of causal directions for pairs of genome sequences with variable lengths. A total of 16,619 high quality complete nucleotide sequences obtained from human hosts from 19 countries were obtained -12,556 from the GISAID initiative's EpiCoV database [47] and the remaining from GenBank [48] . Each sequence was encoded numerically with the following mapping: A=1, C=2, G=3, T=4. Sequences with ambiguous nucleotides were not subject to analysis. A complete table of all sequences used, with accession identifiers and acknowledgements, is available in an online GitHub repository at https://github.com/pranaysy/CausalDiscoveryPaper
All experiments were executed in parallel across multiple cores on a workstation with a 1.8GHz 16-core Intel Xeon Silver 4108 CPU and 32GB of memory, running 64-bit Ubuntu 16.04LTS on kernel version 4.4.
With varying strengths of coupling for coupled AR processes, we observe in Fig. 1 that the accuracies of ETC-P, ETC-E and ERGO improve slightly with coupling while that of ORIGO remains constant and that of LZ-P worsens with coupling. Estimates of accuracy for each value of coupling strength are derived from 1,000 trials. The average accuracy obtained for coupling strengths greater 0.3 is highest for ERGO at 72%, followed by ETC-E and ETC-P at 61% and 58% respectively.
We examine the performance of all models as accuracy against the fraction of decisions each model is forced to make, or decision rates. As shown in Fig. 2 , with stronger coupling overall accuracies improve with decision rates for ETC-E, ETC-P and ERGO. If we focus at the top 20% most decisive coupled processes for which the strength of causal evidence was largest, ETC-P shows an accuracy of 80%, followed by ERGO and ETC-E with accuracies of roughly 75% and 70% respectively. The performance of ORIGO worsens slightly with increasing coupling, while that of LZ-P drops to an accuracy of 0% with coupling strengths greater than 0.3.
We further investigate the performance of all models using estimates of AUROC and AUPRC, summarized across varying coupling strengths in Fig. 3 . We observe an increasing trend in AUROC estimates for ETC-P with coupling strengths, indicating that larger values of causal strengths obtained from ETC-P tend to favor the true causal direction and that this relationship is stronger with stronger coupling. While AUROC for ETC-E, ERGO and ORIGO are approximately constant and slightly better than random for coupling strengths greater than 0.3, and that of LZ-P reaches 0. AUPRC estimates reveal that ETC-E, ETC-P and ERGO are able to correctly identify the true direction with similar performance.
We investigate the accuracy of all five models against the fraction of decisions each model is forced to make. In Fig.  4 , we show the accuracy versus the decision rate for the benchmark Tübingen cause-effect pairs. If we look at all the pairs, we observe that ETC-E and LZ-P infer the correct direction roughly 60% and 58% respectively, followed by 54% for ORIGO and ERGO. ETC-P infers the correct direction for only 42% of all pairs. Considering only those pairs for which the strength of causal evidence was relatively large, we observe that over the top 8% most decisive pairs, ETC-E is 78% accurate and 70% accurate for the top 20% pairs. Similar accuracies are observed for LZ-P.
ROC curves for all models shown in Fig. 5 suggest that causal strength is not very reliable at prediction of causal direction and performance is only slightly better than random. Large estimates of causal strengths do not necessarily correspond to the true direction. Estimates of AUROC, summarized in Table 1 , reflect the same. We do observe that of all models, AUROC as well as AUPRC are highest for ETC-E and LZ-P. The performance of these two proposed models compares favorably with the state-of-the-art causal inference frameworks for continuous real-valued data [34, 6, 49, 9] .
We carried out a novel application of the proposed models for causal inference of the direction of information flow between genomic sequences of the SARS-CoV-2 virus. We evaluated three hypotheses regarding causal discovery in genomic sequences using the ETC-E, ETC-P and LZ-P models.
Firstly, we assessed whether the SARS-CoV-2 consensus sequence (RS) 'causes' all other sequences. 9 countries of 16 for ETC-E, 18 countries of 19 for ETC-P and 11 countries of 16 for LZ-P had at least 5% sequences which admitted this hypothesized direction. Of the three models, higher proportions across countries were generally observed under the ETC-P model.
Secondly, we assessed whether the first SARS-CoV-2 genome sequence isolated in each country (CW) 'causes' all other sequences isolated in that country. 18 countries of 19 for ETC-E, 17 countries of 19 for ETC-P as well as LZ-P had at least 5% sequences which admitted this hypothesized direction. Of the three models, higher proportions across countries were again generally observed under the ETC-P model.
Lastly, we assessed whether the proportions of sequences per country 'caused' by CW were greater than RS. This was the case with 14 countries for ETC-E, 11 for ETC-P and 12 for LZ-P out of 19 countries. Proportions per country obtained under each model are presented in Fig. 6 and estimates provided in Supplementary Table 1 .
Our experiments demonstrate that ETC-P and ETC-E are able to infer causal direction adequately in practice. On synthetic data, both models perform well for moderate to high coupling strengths, while on benchmark data they perform at least as well as the compared models. We observed this performance despite the fact that both these experiments involved information loss due to discretization of continuous real-valued data. We attribute this performance to effective inference of unique context-free grammars (CFGs) by ETC from input symbolic sequences. The inferred grammars likely robustly capture information content of sequences, permitting an effective assessment of information shared and exchanged between sequences by both ETC-P and ETC-E.
The cumulative entropy-based ERGO performs better than other models for synthetic data on the grounds of overall accuracy. Accuracy over top 20% most decisive coupled pairs as well as AUPRC estimates highlight lack of differences in performance between ERGO, ETC-E and ETC-P. AUROC estimates on the other hand reveal an inadequate performance of ERGO, particularly for moderate to strong degrees of coupling. Differences can potentially arise due to the robustness of CCMs over entropy as complexity measures especially for short noisy time series, as well as due to the impact of conditioning. Under the two-part definition of Kolmogorov complexity used for ERGO, the cost of conditioning one sequence using the compressed representation of another sequence is always non-negative. While in the case of ETC-P, conditioning a sequence on the grammar of another sequence may result in either greater or poorer compression (CCM L (y|G x ) can be greater than, less than or equal to CCM L (y)), resulting in positive or negative or zero penalties respectively, since P x→y = CCM L (y|G x ) -CCM L (y). This likely allows ETC-P to capture different patterns of information exchange between sequences.
The systematically poor performance of LZ-P in synthetic experiments suggests an alternate interpretation of the proposed model itself. We argue that due to the algorithmic differences between LZ and ETC, the estimated conditional CCMs (CCM L (y|G x ) and CCM L (x|G y )) may be capturing different information structures. In the ETC-P case, the term CCM L (y|G x ) decomposes into ETC(y|G x ) and ETC(y residual ), where the former term corresponds to the cross-compression or conditioning step. This step halts when the inferred grammar G x can not compress y further. However in the case of LZ-P, CCM L (y|G x ) is estimated as LZ(xy), where xy is the sequence obtained by concatenation of x and y. Here, CCM L (y|G x ) implicitly includes the CCM of x, leading to an interpretation in terms of information accumulation. LZ-P might capture the information accumulated firstly in x since it is parsed before y and then in y but conditioned on x, beyond information contained in y itself (LZ(y)) without any conditioning. The inference of direction now corresponds to greater accumulation of information, leading to results that are aligned with the ground truth of the synthetic experiments. We argue that either the "lesser penalty" or the "higher accumulated information" interpretation may hold true for inference of direction given our proposed framework, and further experiments are needed to investigate and identify conditions suitable for either.
Our unique application of the proposed models for causal inference in SARS-CoV-2 genome sequences offers insights into directional information exchange between sequences. While each of the three models employed provide different causal perspectives, they all relied on inferred generative grammars underlying nucleotide sequences which uniquely characterize them. The results for causal inference are encouraging and showcase an application to genome sequences with very different lengths which does not require sequence alignment and examines sequences globally by decomposing their grammars. We believe these properties to be highly desirable for various applications in bioinformatics and systems biology such as estimating network graphs of sequence or gene interactions through pairwise causal inference, phylogenetics, motif discovery, quantifying dependence and interactions between sequences, among others [50, 51, 52, 42] . These applications may be beneficial to epidemiological studies that rely on temporal information such as contact-tracing, founder effect, epidemic monitoring, evolution of virulence and pathogenicity, etc [40, 53] .
While our experiments spanned synthetic, benchmark and real world data, a more extensive evaluation of performance of the proposed models is needed. In addition to linear AR coupling, linearly as well as non-linearly coupled dynamical systems, such as skew-tent maps, may characterize model performance with greater detail. Robustness to various kinds of noise and scaling, as well as applications to various publicly available benchmark data sets may help elucidate differential robustness, if any, of each of the models. The ETC-E model could particularly benefit from investigations into a more suitable normalization factor, λ, which may yield a better and more consistent performance for causal inference. The overall performance of all models tested is ultimately impacted by the choice of a discretization strategy, and while this remains an open problem, the usage of different methods such as K-Nearest Neighbours (KNNs) [54] , Symbolic Aggregate approXimation (SAX) [55] , etc, may be considered definite improvements.
We presented a new grammar-based information theoretic framework for univariate causal inference using Compression-Complexity Measures (CCMs) and proposed two distinct approaches based on this framework. The penalty and efficacy approaches estimate the cost and effectiveness, respectively, of compressing a sequence using a grammar other than its own. We use these estimates in both directions for a pair of sequences to identify a causal direction based on less penalty and greater efficacy respectively.
For application in practice, we implemented these models using the lossless compressors ETC and LZ which serve as CCMs. The three models -ETC-P, ETC-E and LZ-P allow for reliable causal inference without making any assumptions about the data or its distributions. Empirical evaluation showed that the models are very accurate and perform competitively when compared to the state-of-the-art methods. Further, the fact that these models are essentially parameter-free makes their application in practice very easy compared to methods like Transfer Entropy [4] , or CCC [5] which require careful parameter selection.
We presented a novel application of these models to causal analysis of SARS-CoV-2 genome sequences and found encouraging results which we believe can offer further insights into the disease transmission and sequence evolution. Our models are capable of discovering causal interactions directly from sequences, presenting opportunities for targeting key issues in bioinformatics, genomics and systems biology. As future work, we intend to refine the models presented, evaluate their performance on different causal interactions, and extend them for inference of causal networks.