# Simulation-based Benchmarking for Causal Structure Learning in Gene
Perturbation Experiments

Luka Kovačević
MRC Biostatistics Unit
University of Cambridge
Cambridge, UK
<EMAIL_ADDRESS>
&Izzy Newsham
MRC Biostatistics Unit
University of Cambridge
Cambridge, UK
Sach Mukherjee
DZNE and University of Bonn, Bonn, Germany,
MRC Biostatistics Unit
University of Cambridge, Cambridge, UK
&John Whittaker
MRC Biostatistics Unit
University of Cambridge
Cambridge, UK

###### Abstract

Causal structure learning (CSL) refers to the task of learning causal
relationships from data. Advances in CSL now allow learning of causal graphs
in diverse application domains, which has the potential to facilitate data-
driven causal decision-making. Real-world CSL performance depends on a number
of context-specific factors, including context-specific data distributions and
non-linear dependencies, that are important in practical use-cases. However,
our understanding of how to assess and select CSL methods in specific contexts
remains limited. To address this gap, we present CausalRegNet, a
multiplicative effect structural causal model that allows for generating
observational and interventional data incorporating context-specific
properties, with a focus on the setting of gene perturbation experiments.
Using real-world gene perturbation data, we show that CausalRegNet generates
accurate distributions and scales far better than current simulation
frameworks. We illustrate the use of CausalRegNet in assessing CSL methods in
the context of interventional experiments in biology.

###  1 Introduction

The aim of causal structure learning (CSL) is to infer causal relationships
between variables, typically in the form of a directed graph. CSL has been
extensively studied and a wide array of methods are available in the
literature [1]. Recent years have seen considerable theoretical and
methodological advances in CSL, including reformulation via continuous
optimization [2], neural methods [3, 4] and emerging, end-to-end supervised
deep learning-based CSL approaches [5, 6].

However, despite these exciting developments, CSL in scientific applications
remains challenging due to a number of factors, including problem dimension,
noise, data limitations and latent variables. Although many CSL methods are
supported by theory, in practice it is challenging to verify underlying
assumptions. Thus, it is often difficult to assess the likely efficacy of a
CSL approach in a real problem context. This means that in many contexts –
often involving arduous and costly experiments – it remains nontrivial to
judge whether a given CSL approach is well suited to the setting.

The application of CSL to problems in biomedicine manifests many such
problems. Advances in cell and molecular biology have dramatically improved
our ability to perform interventional experiments at scale. Notably, the
combination of gene editing technologies (e.g. CRISPR) with sequencing
protocols [7] has enabled measurement of tens of thousands of variables under
thousands of interventions, with recent experiments spanning millions of cells
[8]. These developments are in the process of transforming biological and
biomedical discovery, enabling causal experiments at an unprecedented scale.
Large-scale gene knockout experiments in particular, promise to inform drug
target identification by allowing for inference of the effects of individual
and combinatorial genetic interventions. An emerging body of work has focused
on CSL for biological data [among others 9, 4, 6], but for a given problem
setting it remains practically challenging to understand efficacy and
therefore to establish reliable CSL-based workflows.

###### Related Work.

The need to better understand CSL performance characteristics has inspired a
line of work in benchmarking methods [10, 11]. While a number of papers
include evaluation using real biological data [9, 4, 6], these results are
inevitably specific to the particular real datasets used and relative
performance of CSL methods can vary significantly depending on context. This
has motivated work on simulation engines of various kinds. In the biological
domain, Single-cell Expression of Genes In-silico (SERGIO) [12] was developed
as an approach based on stochastic differential equations. More recently,
GRouNdGAN [13] was proposed as a technique for generating single-cell RNA
sequencing (scRNA-seq) data from a given causal graph. CausalBench [11] has
also been proposed as a framework for evaluating structure learners directly
from real-world data using gene regulatory networks derived from literature
and statistical hypothesis testing.

Large-scale data simulation is an important part of the deep learning CSL
method reported in [5]. Nonetheless, our approach could be used to generate
data to train, or pre-train, supervised CSL methods, a point we return to in
the Discussion.

###### Contributions.

We propose CausalRegNet, an approach for generating realistic synthetic data
from an underlying directed acyclic graph (DAG). CausalRegNet allows
generation of both observational and interventional data whilst respecting
user-specified statistical properties and conforming to a number of basic
desiderata. Crucially, CausalRegNet is highly scalable, allowing study of
problems of the size now seen in contemporary experiments and data acquisition
pipelines. This simulator can be used to study the performance of CSL methods
in a context-specific manner, by fixing key aspects of the problem setting.

Our work has several novel contributions:

  1. (i)

We present a novel tool called CausalRegNet111Code is available at
https://github.com/luka-kovacevic/causalregnet. for simulating realistic
scRNA-seq data with respect to a causal graph;

  2. (ii)

We validate CausalRegnet using data from real world large-scale interventional
experiments in biology;

  3. (iii)

We demonstrate how CausalRegNet can be used to benchmark CSL algorithms in
practice.

The methods, software and workflow we provide allow for assessment of CSL
methods in the context of real-world problems. Our aim is not to provide a
comprehensive empirical benchmarking study but rather demonstrate how such a
benchmark can be carefully devised whilst taking into account many context-
specific factors and criteria. Our simulator generates data that behaves
similarly to data from (specific) real experiments, allowing for inferences
that could not be made using an entirely generic simulation model.

###  2 Background

In Section 2.1 we introduce the general background knowledge needed for
structural causal models (SCMs). This is followed by the introduction of
additive noise models (ANMs) in Section 2.2, which are a specific subclass of
SCMs that are commonly used to generate data for comparative benchmarks of CSL
methods.

####  2.1 Structural Causal Models

######  Definition 2.1 (Structural Causal Model; [14]).

A structural causal model (SCM) $\mathfrak{C}=(S,\mathcal{P}_{N})$ is a
collection $S$ of $d$ structural assignments corresponding to each
$X_{j}\in\mathbf{X}=\\{X_{1},\ldots,X_{d}\\}$,

$X_{j}:=f_{j}(X_{\text{pa}(j)},N_{j}),\;\;\;\;\;\;j=1,\ldots,d,$ (1)

where $X_{\text{pa}(j)}\subseteq\mathbf{X}\backslash\\{X_{j}\\}$ are parents
of $X_{j}$, $\mathcal{P}_{N}$ is a joint distribution over noise variables,
which we require to be jointly independent, such that we can write
$\mathcal{P}_{N}$ as a product of the marginal densities of $N_{j}$ for
$j=1,\ldots,N$. If $X_{i}\in X_{\text{pa}(j)}$ then $X_{j}$ is a child of
$X_{i}$ and $X_{i}$ is a parent of $X_{j}$.

The causal graph $\mathcal{G}=(V,E)$ is composed of vertices $V=[d]$ and edges
$E=\\{i\rightarrow j:X_{i}\in X_{\text{pa}(j)}\\}$. We assume that
$\mathcal{G}$ has no paths of the form $j\rightarrow\ldots\rightarrow i$, that
is, $\mathcal{G}$ is acyclic and there are no paths along directed edges that
lead from a child to one of its parents. For any node $j$ if there exists a
(directed) path $i\rightarrow\ldots\rightarrow j$ along directed edges then
node $i$ is an ancestor of $j$. Conversely, node $j$ is a descendant of node
$i$.

Thus, the observational distribution over $\mathbf{X}$ implied by an SCM is
such that each $X_{j}=f_{j}(X_{\text{pa}(j)},N_{j})$ in distribution. Assuming
there are no interventions, the distribution entailed by the SCM is the
observational distribution, also more generally called an entailed
distribution.

###### Interventions.

In the SCM framework, interventions can be performed on individual nodes or
groups of nodes, where the set of intervention targets is denoted by
$I\subseteq V$. When considering more than one intervention $I$, the
interventional family $\mathcal{I}:=\\{I_{1},\ldots,I_{K}\\}$ contains $K$
distinct intervention settings. The observational setting can be denoted as
$I^{\emptyset}$. The distribution entailed by the SCM following an
intervention is referred to as an interventional distribution.

Interventions are difficult to carry out in practice and at best limited in
scope in many settings.

However, in many molecular biology settings interventional experiments are
possible, e.g. by interventions targeted at specific molecular mechanisms.
Ongoing technological developments, including methods rooted in CRISPR and
related gene editing protocols, are now enabling large-scale interventional
experiments including those further described in Appendix A.

####  2.2 Additive Noise Models

######  Definition 2.2 (Additive Noise Model; [15]).

An SCM with structural assignments of the form,

$X_{j}:=f_{j}(X_{\text{pa}(j)})+N_{j}$ (2)

where $f_{j}(\cdot)$ is an arbitrary function and the noise variables $N_{j}$
are jointly independent, is an additive noise model (ANM).

The case where the function $f_{j}(\cdot)$ is linear and the noise variables
are distributed as $N_{j}\sim\mathcal{N}(0,\sigma^{2}_{j})$ yields the general
family of linear Gaussian ANMs. We refer to these simply as Gaussian ANMs
throughout.

Synthetic data generated through linear ANMs has previously been shown to be
vulnerable to CSL methods that exploit varsortability to improve causal
structure learning performance.

######  Definition 2.3 (Varsortability; [16]).

Varsortability, denoted by $v$, is a measure of the agreement between the
order of increasing marginal variance and the causal order,

$v:=\frac{\sum_{k=1}^{d-1}\sum_{i\rightarrow j\in
E^{k}}\text{\frenchspacing\mbox{increasing}}(\text{Var}(X_{i}),\text{Var}(X_{j}))}{\sum_{k=1}^{d-1}\sum_{i\rightarrow
j\in E^{k}}1},$ (3) $\text{where
\frenchspacing\mbox{increasing}}(a,b)=\begin{cases}1&a<b,\\\ 1/2&a=b,\\\
0&a>b.\end{cases}$ (4)

Methods that score close to 1 or close to 0 are varsortable, meaning the
correct causal ordering can be obtained by ordering the nodes by ascending or
descending variance, respectively.

###  3 Methodology

To translate recent developments in CSL to the real-world inference of causal
structures, it is essential that we understand exactly when CSL algorithms
fail and why. This type of understanding requires both benchmarking studies
and theoretical work that accounts for context-specific data properties.
Motivated by the foregoing, we focus on context-specific simulation for CSL.
In Section 3.1, we begin by deriving simulator desiderata that are important
to a context-specific simulation tool for large-scale gene perturbation
experiments. Section 3.2 describes the model developed based on desiderata,
including an example of conditions relating to the functional relationship
between parents and children that ensures our model adheres to domain
knowledge of observational and interventional relationships present in scRNA-
seq data.

####  3.1 Desiderata

The desiderata applied to develop CausalRegNet are as follows:

  1. I.

Interpretability: Simulators that are parametrically interpretable are easier
to use and thus allow for simpler reasoning when it comes to comparing models.

  2. II.

Scalability: Contemporary datasets can be very high dimensional, motivating a
need for scalability. Scalability is particularly important for simulation
datasets that are generated from a large set of graph configurations as needed
for model validation or supervised model training.

  3. III.

Low varsortability: Data generated by simple ANMs is highly varsortable.
Previous work by [16] has established that several CSL algorithms are able to
exploit varsortability to yield good benchmark results that are driven by this
property. Thus far, there is no evidence that real-world causal structures
lead to varsortable data, and hence, our simulator of choice should reflect
this fact.

  4. IV.

Distributional Properties: Agreement between the properties of empirical and
synthetic distributions is crucial to ensuring that novel algorithms perform
as well in the real world as they do in theory.

These desiderata guide the design decisions described in Section 3.2 where a
different set of desiderata may have led to a different model.

####  3.2 CausalRegNet

We propose a multiplicative effect SCM that accurately reflects the
interaction dynamics of single-cell gene expression. Once parameterised,
samples can be drawn efficiently from this simulator.

#####  3.2.1 Node-wise Distribution

For each node $j$, the random variable associated with this node is
distributed according to,

$X_{j}\sim\text{NegBin}(\mu_{j},\sigma_{j}^{2}),$ (5)

where $\mu_{j}$ is the mean expression and $\sigma^{2}_{j}$ is the variance.
The parameters are defined as

$\mu_{j}=\mu_{j}^{0}\cdot
f_{j}^{\text{reg}}(X_{\text{pa}(j)};\Theta_{j}),\;\;\;\;\text{ and
}\;\;\;\;\sigma_{j}^{2}=\mu_{j}(1+\mu_{j}/\theta_{j}),$ (6)

where $\mu_{j}^{0}$ is the observational mean and $f_{j}^{\text{reg}}(\cdot)$
represents the regulatory effect of the parents of $X_{j}$ on $X_{j}$. A mean-
variance relationship is imposed via the variance, $\sigma_{j}^{2}$, which is
determined by the mean, $\mu_{j}$, and the inverse dispersion parameter,
$\theta_{j}$. This functional form is based on well-established domain
knowledge motivated by generic and well-known properties of biochemical
interplay [17].

#####  3.2.2 Regulatory Function

The regulatory effect function takes a sigmoidal form,

$f_{j}^{\text{reg}}(X_{\text{pa}(j)};\Theta_{j})=\frac{\alpha_{j}}{1+\exp{\left\\{-\gamma_{j}(w(x_{\text{pa}(j)})+b_{j})\right\\}}},$
(7)

where $\alpha_{j}$ is the maximum regulatory effect of $X_{\text{pa}(j)}$ on
$X_{i}$, and $w(x_{\text{pa}(j)})$ is any function that aggregates the
regulatory effect of parents. This regulatory function introduces non-
linearity to the relationship between each node and its parents.

For simplicity, we use $w(x_{\text{pa}(j)})=\sum_{i\in\text{pa}(j)}w_{ij}\cdot
x_{i}/\mu_{i}^{0}$, where $W\in\mathbb{R}^{n\times n}$ is a weighted adjacency
matrix that defines the strength of relationship between $X_{i}$ and $X_{j}$
and $[W]_{ij}=w_{ij}$ is the edge weight for the edge going from $X_{i}$ to
$X_{j}$.

#####  3.2.3 Model Specification

The parameters $\gamma_{j}$ and $b_{j}$ are calibrated according to the
aggregation function $w(\cdot)$, maximum regulatory effect $\alpha_{j}>1$, and
a minimum regulatory effect $0<\beta_{j}<1$ parameter. To do this, we define
the following value conditions:

  1. I

Baseline Expression Condition (BEC):

$\text{if }\forall i\in\text{pa}(j),X_{i}=0,\text{ then
}f_{j}^{\text{reg}}(X_{\text{pa}(j)})=\beta_{j},$

  2. II

Observational Expression Condition (OEC):

$\text{if }\forall i\in\text{pa}(j),X_{i}=\mu_{i}^{0}\text{ then
}f_{j}^{\text{reg}}(X_{\text{pa}(j)})=1.$

In other words, we first have that when the expression of each of the parents
of node $j$ are at zero, the regulatory effect will be at the baseline
$\beta_{j}$. Conversely, when each of the parents of $j$ are expressed at
their observational mean, then the regulatory effect is equal to 1, implying
$\mu_{j}=\mu_{j}^{0}\cdot 1=\mu_{j}^{0}$. That is, $X_{j}$ is at its
observational mean.

These value conditions allow us to freely specify the aggregation function
$w(\cdot)$ and calibrate the regulatory effect function to guarantee maximum
and minimum effects. Furthermore, these conditions are sufficient to yield
closed-form solutions for $\gamma_{j}$ and $b_{j}$ given the parameters
$\\{\alpha_{j},\beta_{j},w(\cdot)\\}$.

Thus, given that the node-wise distribution, regulatory effect function and
model specification method have clear biological interpretations, CausalRegNet
satisfies Desiderata I.

#####  3.2.4 Example: Model Specification for Linear Aggregation with Mean-
normalisation

Assuming the linear aggregation function introduced in Section 3.2.2, we
demonstrate how value conditions can be used to derive the parameters
$\gamma_{j}$ and $b_{j}$. Applying OEC and BEC from Section 3.2.3, we have
that $w(x_{\text{pa}(j)})=0$ when $x_{i}=0$, $\forall i\in\text{pa(j)}$ and
that $w(x_{\text{pa}(j)})=\sum_{i\in\text{pa}(j)}w_{ij}$ when
$x_{i}=\mu_{i}^{0}$, $\forall i\in\text{pa}(j)$, respectively. Solving, this
yields,

$b_{j}=-\frac{1}{\gamma_{j}}\left(\log{\left(\alpha_{j}-1\right)}+\gamma_{j}w^{\prime}_{j}\right)\;\;\;\;\text{and,}$
(8)
$\gamma_{j}=\frac{1}{w^{\prime}_{j}}\left[\log{\left(\frac{\alpha_{j}}{\beta_{j}}-1\right)}-\log{\left(\alpha_{j}-1\right)}\right],$
(9)

where $w^{\prime}_{j}=\sum_{i\in\text{pa}(j)}w_{ij}$ and $w^{\prime}_{j}\neq
0$. Assuming that $w^{\prime}_{j}=1$ and setting the remaining parameters as
$\\{\alpha_{j}=2,\beta_{j}=0.1\\}$ the calibrated regulatory effect function
is shown in Figure 1. The points in this figure highlight the OEC and BEC
being satisfied on the regulatory effect curve.

Figure 1: Example calibrated regulatory effect function. The red line
indicates the cut-off, the regulatory effect will never go left of this line.

#####  3.2.5 Adjusted Regulatory Function

When simulating nodes with excessively small observational means
$\mu^{0}_{j}$, standard calibration as shown in Figure 1 can lead to overly
large aggregated regulatory effects. For example, suppose that we want to
simulate a child of gene MTOR from [8] with $\hat{\mu}^{0}_{j}\approx 0.25$
and $\hat{\theta}_{j}\approx 3.38$ with $f^{\text{reg}}$ shown in Figure 1.
The distribution of the expression of MTOR is shown in Figure 5(a), where it
varies between 0 and 4. With such expression rates, we have that
$x_{i}=0,1,2,3,4$, which corresponds to $x_{i}/\mu_{i}^{0}=0,4,8,12,16$,
respectively. This regulatory effect is much larger than the range of values
for which the function in Figure 1 varies. That is, for any $x_{i}>1$,
$f^{\text{reg}}_{j}(X_{\text{pa}(j)})\approx 2$. Thus, we propose an adjusted
regulatory function with an additional regulatory adjustment constant $a_{j}$,

$f^{\text{reg}}_{j}(X_{\text{pa}(j)};\Theta_{j})=\frac{\alpha_{j}}{1+\exp{\\{-\gamma_{j}(w_{\text{adj}}(x_{\text{pa}(j)})+b_{\text{adj}})\\}}}$
(10)

where

$b_{\text{adj}}=b_{j}-a_{j},\;\;\;\;\text{and}\;\;\;\;w_{\text{adj}}(x_{\text{pa}(j)})=w(x_{\text{pa}(j)})+a_{j}.$
(11)

This additional hyperparameter now allows for full expression of regulatory
effects even when $\mu_{i}^{0}$ is very small. This adjusted regulatory effect
function is used for simulations throughout this work.

#####  3.2.6 Node-wise Distribution Fitting

By parametrising the distribution of each node as a Negative Binomial, we
enable node-wise fitting of simulated distributions to empirical marginal
distributions in real-world data. We fit the observational mean and inverse
dispersion via the Negative Binomial likelihood [17]. The implemented fitting
scheme uses the L-BFGS-B to find a locally optimal solution for each node
independently. Further details can be seen in Appendix C.

###  4 Empirical Results

The development of CausalRegnet was motivated by the lack of agreement between
interventional distributions generated by SERGIO and domain-knowledge in
biology. We explore this motivation in Section 4.1 alongside results from
CausalRegNet. Guided by our desiderata, we then compare the scalability and
varsortability of CausalRegNet, SERGIO and the Gaussian ANM in Sections
4.2-4.3, respectively. Finally, in Section 4.4 we examine the agreement
between CausalRegNet and real-world perturb-seq data. The full parameter
settings used for simulation in this section are given in Appendix D.2.

####  4.1 Motivation: Interventional Distributions under Changing Structures

To understand simulator behaviour under intervention, we estimate the average
treatment effects (ATEs) of the intervention $I=\\{X_{0}\\}$, on $X_{1}$,
denoted

$ATE(X_{0},X_{1})=\mathbb{E}\left[X_{1}^{\emptyset}\right]-\mathbb{E}\left[X_{1}^{\textit{do}(X_{0}=0)}\right],$
(12)

with varying causal structures, where $\emptyset$ denotes the observational
setting.

${X_{0}}$$c_{1}$$X_{1}$$c_{k}$… (a) More parents

${X_{0}}$$c_{1}$$c_{k}$$X_{1}$… (b) More mediators

Figure 2: We examine the treatment effect of $X_{0}$ on $X_{1}$ as (a) the
number of parents of $X_{1}$ increases and (b) as the number of mediators
between $X_{0}$ and $X_{1}$ increases (i.e. for $k=0,\ldots,9$ under (a) and
(b)).

First, we study the ATE as the number of parents of a node increases. That is,
we are interested in how the distribution of a random variable changes as the
number of parents increases. Figure 2(a) illustrates this setting, where the
number of parents, $k$, varies.

In Figure 3a, simulating with SERGIO gives an average ATE of 2.5 across
simulations, whilst the mean expression of $X_{1}$ is rising. This suggests
that the constant ATE is facilitated by the increasing expression level,
allowing each additional parent to have an ATE of 2.5. Conversely,
CausalRegNet with the mean-normalised aggregation function applied in Section
3.2.3 ensures that each regulator has equal contribution, and thus, as the
number of parents increases, the observed ATE tends towards zero. Although it
is not clear what aggregation function is most faithful to the underlying
biology, the modularity of CausalRegNet means that the aggregation function
can be adapted to emerging domain knowledge. For both CausalRegNet and SERGIO,
the ATE is partially determined by the weighted adjacency matrix, such that
parents with greater edge weights will determine more of the regulatory effect
and will result in different ATEs. This can be directly altered when using
either simulator.

Figure 3: (a) shows absolute ATE with increasing numbers of parents. In (b),
the mean observational expression of $X_{1}$ against number of parents. (c)
shows the absolute ATE falls for both simulators as the number of mediators
increases.

Next, we consider the case where the number of nodes mediating the effect of
$X_{0}$ on $X_{1}$ increases. Figure 2(b) illustrates this setting, where
$k=0,\ldots,9$. Figure 3c shows that as the number of mediators increases, the
ATE falls. This is intuitively expected since ancestors further away may
plausibly lead to smaller changes in downstream nodes. Thus, both simulators
behave as expected under this case.

####  4.2 Scalability

Figure 4a shows the time taken to simulate data from predetermined graphs of
size 3 to 10,000. The scalability of CausalRegNet is shown to be comparable to
the simple Gaussian ANM, which can generate data from graphs of size 1,000
within less than 10 seconds. For the same graphs, SERGIO takes upwards of
10,000 seconds or 2.78 hours. Moreover, SERGIO takes over 100,000 seconds or
24 hours to simulate a causal graph of 5,000 or more nodes, which makes it
practically very difficult to simulate such large datasets.

The poor scalability of SERGIO is caused by the need to reach a steady state
before samples can be drawn. Methods based on SCMs such as CausalRegNet or any
ANM are defined by parametric probability distributions and can thus be be
drawn from instantaneously. The simulation time in this instance is purely
associated with the time needed to traverse the causal graph according to the
topological ordering of nodes.

Thus, CausalRegNet satisfies Desiderata II as it is scalable to the dimension
of real-world experiments.

Figure 4: (a) Simulation time for data from graphs of increasing size. Data is
generated from the same causal graph structures across methods. Computational
time for a single simulation with SERGIO with 5,000 nodes or more exceeds 24
hours and so is not included. Mean varsortability of data generated from (b) a
causal chain and (c) a causal graph structure by each simulator with 95%
confidence intervals.

####  4.3 Varsortability

To examine the varsortability of data generated by each simulator we generated
data from causal chain and causal graph structures. A causal chain is a causal
graph such that only one directed path of the form
$X_{1}\rightarrow\ldots\rightarrow X_{n}$ exists, starting from the first node
in the chain, $X_{1}$, to the last $X_{n}$. For CausalRegNet, we use nodes
fitted to real gene expression data from [8] by randomly selecting genes from
the set obtained in Appendix D.4. Parametric settings for simulating from
SERGIO can be found in Appendix D.2. Parameters used to generate causal graphs
used by both CausalRegNet and SERGIO are in Appendix D.3.1.

The results for the causal chain simulations are shown in Figure 4b. For
causal chains of length 10 or more, CausalRegNet and SERGIO generate near
perfectly non-varsortable data. In Figure 4c CausalRegNet generates data with
very low varsortability and outperforms both SERGIO and the Gaussian-ANM in
this respect. Further, Table 1 shows the average varsortability of data
generated by the causal chains and graphs with the greatest number of nodes
(i.e. $d=\\{30,50\\}$, respectively), with CausalRegNet achieving
significantly lower varsortability than both SERGIO and the Gaussian ANM.

CausalRegNet therefore satisfies Desiderata III as it generates data with low
varsortability.

Table 1: Varsortability of data generated by the largest causal chains ($d=30$) and causal graphs ($d=50$) with 95% confidence intervals. Simulator | Causal Chain Varsort. | Causal Graph Varsort.
---|---|---
Gaussian ANM | $0.863\pm 0.076$ | $0.986\pm 0.005$
SERGIO | $\mathbf{0.488\pm 0.034}$ | $0.873\pm 0.019$
CausalRegNet | $\mathbf{0.477\pm 0.008}$ | $\mathbf{0.492\pm 0.008}$

####  4.4 Comparison to Experimental Data

To conduct an empirical comparison between the data generated by CausalRegNet
and real-world data, we collect 100 genes related to cancer pathways (see
Appendix D.4) for the most commonly referenced genes that were both measured
and perturbed in the experimental dataset from [8].

Given this set of 100 genes, we aim to simulate graphs of size
$d=\\{3,10,100\\}$ with CausalRegNet to assess whether the simulated
distributions are similar to those observed in [8]. Using the distribution
fitting scheme outlined in Section 3.2.6, each node in the DAG is fitted to
simulate the expression of a gene from the original set of 100 genes.

Figure 5a shows the real and simulated nodes for 3 genes selected from [8]. We
measure the distance between simulated and real distributions using the
Wasserstein distance (WD) [18, 19]. The overlap between simulated and real
distributions is $<0.2$ for all nodes. Figure 5b shows the WDs for graphs of
size $d=100$. As a comparative baseline, we shuffle the gene labels such that
each synthetic distribution is compared to the expression of a random gene.
From Figure 5b, clearly the comparison to the true distribution leads to lower
WDs. For both true and random label assignments, there are also outliers in
the WDs, shown in Figure 8 (Appendix D.6), however, these account for a small
number of nodes. To assess the synthetic interventional distributions, since
the true causal graph is unknown we cannot compare distributions directly.
Instead, in Figures 5c-d we compare the distribution of interventional effects
between the real and simulated data. We can see that we obtain good coverage
with some excessive $0$ effects in 5c, however, modifying $\alpha_{j}$ leads
to a smaller number of $0$ as shown in Figure 5d.

Figure 5: (a) Comparison between real and simulated data for nodes fitted to
genes in [8] in DAG with 3 nodes. (b) In green, WD for simulated marginal
distributions fitted to [8] compared to the true distributions. In purple, the
distribution labels are shuffled such that each simulated marginal is compared
to a random empirical marginal, giving a baseline for comparison. The
distribution of interventional effects with (c) $\alpha_{j}=2$ and (d)
$\alpha_{j}=5$ for each node.

Hence, CausalRegNet also satisfies Desiderata IV as it generates observational
and interventional distributions that are comparable to real-world data.

####  4.5 Benchmark Study

The goal of this work is to introduce a context-specific simlulator for use in
studying CSL. To demonstrate its use, we carried out an illustrative study
benchmarking some CSL algorithms commonly used as baselines (see Appendix
D.3). Here, we see that in the specific context used, the various methods
perform poorly. We emphasise that fully leveraging our simulators for
benchmarking per se goes beyond the scope of this paper and would be a key
direction for future work.

###  5 Discussion

In this work, we proposed CausalRegNet, a simple yet flexible SCM-based
approach to simulation in the context of large-scale gene perturbation
screens. We found that CausalRegNet led to significantly improved scalability
over existing approaches whilst coping with the issue of varsortability and
providing synthetic data matched to key properties of real data distributions.
Recent work [20] has shown that varsortability corresponds to a necessary
assumption for some ANMs, as shown by [21], thus CausalRegNet again provides
an approach for comparing methods with varying assumptions under a fair,
realistic setting. In addition to its use for benchmarking and context-
specific assessment, we think CausalRegNet can be a useful tool in training
machine learning approaches for causal discovery and inference by allowing for
realistic and rapid simulation. Rapid simulation frameworks can also
contribute to the design of real-world experiments, including active learning
schemes [22, 23, 24].

Individual existing simulators satisfy some but not all of the desiderata set
out in this work. SERGIO is able to generate accurate observational
distributions but lacks the ability to realistically simulate interventions
and scalability to large problems. GRouNdGAN generates accurate observational
distributions and is scalable but lacks interpretability and has not yet been
assessed on its ability to simulate interventions. GRouNdGAN also requires
training and validation on a reference dataset to accurately simulate data,
which can be computationally expensive. In contrast, we provide a library of
parameters derived from [8] that allows simulations to be carried out
parsimoniously.

###### Limitations & Future Work.

Several properties of the framework and tool here require further work to
facilitate development of effective context-specific benchmarking tools. The
simplicity of the SCM-based approach means that temporal dynamics cannot be
simulated, which is possible in less scalable tools such as SERGIO and
approaches that allow time-dependent dynamics simulation will be important in
some applications of CSL. Further work is also needed to explore alternative
parameterisations of the aggregation functions $w(\cdot)$. A key direction for
future empirical work will be to conduct a comprehensive benchmarking study
using CausalRegNet to compare the performance of a library of existing CSL
algorithms.

### Acknowledgments and Disclosure of Funding

LK acknowledges the support of an Ivan D Jankovic Studentship at Clare Hall,
University of Cambridge. This work was funded by UKRI Programme Grant
MC_UU_00002/18. For the purpose of open access, the author has applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising.

### References

  * [1] Christina Heinze-Deml, Marloes H Maathuis, and Nicolai Meinshausen.  Causal structure learning.  Annual Review of Statistics and Its Application, 5:371–391, 2018.
  * [2] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing.  Dags with no tears: Continuous optimization for structure learning.  Advances in neural information processing systems, 31, 2018.
  * [3] Phillip Lippe, Taco Cohen, and Efstratios Gavves.  Efficient neural causal discovery without acyclicity constraints.  arXiv preprint arXiv:2107.10483, 2021.
  * [4] Romain Lopez, Jan-Christian Hütter, Jonathan Pritchard, and Aviv Regev.  Large-scale differentiable causal discovery of factor graphs.  Advances in Neural Information Processing Systems, 35:19290–19303, 2022.
  * [5] Nan Rosemary Ke, Silvia Chiappa, Jane Wang, Anirudh Goyal, Jorg Bornschein, Melanie Rey, Theophane Weber, Matthew Botvinic, Michael Mozer, and Danilo Jimenez Rezende.  Learning to induce causal structure.  arXiv preprint arXiv:2204.04875, 2022.
  * [6] Kai Lagemann, Christian Lagemann, Bernd Taschler, and Sach Mukherjee.  Deep learning of causal structures in high dimensions under data limitations.  Nature Machine Intelligence, 5(11):1306–1316, 2023.
  * [7] Atray Dixit, Oren Parnas, Biyu Li, Jenny Chen, Charles P Fulco, Livnat Jerby-Arnon, Nemanja D Marjanovic, Danielle Dionne, Tyler Burks, Raktima Raychowdhury, et al.  Perturb-seq: dissecting molecular circuits with scalable single-cell rna profiling of pooled genetic screens.  cell, 167(7):1853–1866, 2016.
  * [8] Joseph M Replogle, Reuben A Saunders, Angela N Pogson, Jeffrey A Hussmann, Alexander Lenail, Alina Guna, Lauren Mascibroda, Eric J Wagner, Karen Adelman, Gila Lithwick-Yanai, et al.  Mapping information-rich genotype-phenotype landscapes with genome-scale perturb-seq.  Cell, 185(14):2559–2575, 2022.
  * [9] Steven M Hill, Laura M Heiser, Thomas Cokelaer, Michael Unger, Nicole K Nesser, Daniel E Carlin, Yang Zhang, Artem Sokolov, Evan O Paull, Chris K Wong, et al.  Inferring causal molecular networks: empirical assessment through a community-based effort.  Nature methods, 13(4):310–318, 2016.
  * [10] Marco Eigenmann, Sach Mukherjee, and Marloes Maathuis.  Evaluation of causal structure learning algorithms via risk estimation.  In Conference on Uncertainty in Artificial Intelligence, pages 151–160. PMLR, 2020.
  * [11] Mathieu Chevalley, Yusuf Roohani, Arash Mehrjou, Jure Leskovec, and Patrick Schwab.  Causalbench: A large-scale benchmark for network inference from single-cell perturbation data.  arXiv preprint arXiv:2210.17283, 2022.
  * [12] Payam Dibaeinia and Saurabh Sinha.  Sergio: a single-cell expression simulator guided by gene regulatory networks.  Cell systems, 11(3):252–271, 2020.
  * [13] Yazdan Zinati, Abdulrahman Takiddeen, and Amin Emad.  Groundgan: Grn-guided simulation of single-cell rna-seq data using causal generative adversarial networks.  Nature Communications, 15(1):4055, 2024.
  * [14] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf.  Elements of causal inference: foundations and learning algorithms.  The MIT Press, 2017.
  * [15] Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf.  Nonlinear causal discovery with additive noise models.  Advances in neural information processing systems, 21, 2008.
  * [16] Alexander Reisach, Christof Seiler, and Sebastian Weichwald.  Beware of the simulated dag! causal discovery benchmarks may be easy to game.  Advances in Neural Information Processing Systems, 34:27772–27784, 2021.
  * [17] Davide Risso, Fanny Perraudeau, Svetlana Gribkova, Sandrine Dudoit, and Jean-Philippe Vert.  A general and flexible method for signal extraction from single-cell rna-seq data.  Nature communications, 9(1):284, 2018.
  * [18] Aaditya Ramdas, Nicolás García Trillos, and Marco Cuturi.  On wasserstein two-sample testing and related families of nonparametric tests.  Entropy, 19(2):47, 2017.
  * [19] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272, 2020.
  * [20] Ignavier Ng, Biwei Huang, and Kun Zhang.  Structure learning with continuous optimization: A sober look and beyond.  In Causal Learning and Reasoning, pages 71–105. PMLR, 2024.
  * [21] Po-Ling Loh and Peter Bühlmann.  High-dimensional learning of linear causal networks via inverse covariance estimation.  The Journal of Machine Learning Research, 15(1):3065–3105, 2014.
  * [22] Scott Sussex, Caroline Uhler, and Andreas Krause.  Near-optimal multi-perturbation experimental design for causal structure learning.  Advances in Neural Information Processing Systems, 34:777–788, 2021.
  * [23] Panagiotis Tigas, Yashas Annadani, Andrew Jesson, Bernhard Schölkopf, Yarin Gal, and Stefan Bauer.  Interventions, where and how? experimental design for causal models at scale.  Advances in Neural Information Processing Systems, 35:24130–24143, 2022.
  * [24] Panagiotis Tigas, Yashas Annadani, Desi R Ivanova, Andrew Jesson, Yarin Gal, Adam Foster, and Stefan Bauer.  Differentiable multi-target causal bayesian experimental design.  In International Conference on Machine Learning, pages 34263–34279. PMLR, 2023.
  * [25] Alexander Hägele, Jonas Rothfuss, Lars Lorch, Vignesh Ram Somnath, Bernhard Schölkopf, and Andreas Krause.  Bacadi: Bayesian causal discovery with unknown interventions.  In International Conference on Artificial Intelligence and Statistics, pages 1411–1436. PMLR, 2023.
  * [26] Peter Spirtes, Clark N Glymour, and Richard Scheines.  Causation, prediction, and search.  MIT press, 2000.
  * [27] Liam Solus, Yuhao Wang, and Caroline Uhler.  Consistency guarantees for greedy permutation-based causal inference algorithms.  Biometrika, 108(4):795–814, 2021.
  * [28] Yuhao Wang, Liam Solus, Karren Yang, and Caroline Uhler.  Permutation-based causal inference algorithms with interventions.  Advances in Neural Information Processing Systems, 30, 2017.
  * [29] Jonathan Kans.  Entrez direct: E-utilities on the unix command line.  Entrez Programming Utilities Help [Internet], Apr 2013.  [Updated 2024 Apr 4].
  * [30] Liviu Pirvan and Shamith A Samarajiwa.  Pangaea: A modular and extensible collection of tools for mining context dependent gene relationships from the biomedical literature.  bioRxiv, pages 2020–04, 2020.

## Appendices

In Appendix A, we introduce single-cell RNA sequencing experiments and the
biology necessary for CausalRegNet. Appendix B includes details on SERGIO
further to those included in the main paper. The likelihood used to fit the
negative binomial distributions of CausalRegNet to real data is given in
Appendix C. Finally, Appendix D includes a benchmark on commonly used CSL
algorithms with some additional empirical evaluations of CausalRegNet.

###  Appendix A Introduction to Single-cell RNA Sequencing

Genes are DNA sequences that are transcribed into messenger RNA (mRNA) before
being turned into proteins. The level of mRNA corresponding to each gene
within a particular cell is referred to as the gene expression and is a
measure of the activity of that gene in that specific context. Gene expression
is highly plastic, varying according to the type of cell, tissue, time and
internal and external stimuli, including the expression of other genes. Many
disease associated DNA variants lie outside genes but seem to exert their
effects by influencing RNA expression at genes. Gene expression can be
experimentally measured within individual cells, in single-cell RNA sequencing
(scRNA-seq) experiments, or aggregated across cells within a tissue sample, in
bulk RNA-seq experiments. This in turn allows for the study of the regulatory
structure governing gene expression, enabling us to identify and analyse gene-
gene interactions.

Experimental techniques for collecting scRNA-seq data introduce additional
technical noise, which obscures the already noisy biological signal present in
scRNA-seq data. Existing statistical approaches for studying scRNA-seq data
often assume that data follows a (zero-inflated) Negative Binomial
distribution and have been shown to be effective for signal extraction and
distilling biological and technical noise [17].

###  Appendix B Further details on SERGIO

SERGIO (single-cell expression of genes in-silico) is a technique for
simulating scRNA-seq dynamics using stochastic differential equations (SDEs).
Given a GRN structure, SERGIO simulates gene expression dynamics as a function
of its regulators. It does so using the chemical Langevin equation (CLE)
(cite) where the time-course mRNA concentration of gene $i$ is given by:

$\frac{dx_{i}}{dt}=P_{i}(t)-\lambda_{i}x_{i}(t)+q_{i}\left(\alpha\sqrt{P_{i}(t)}+\beta\sqrt{\lambda_{i}x_{i}(t)}\right)$
(13)

where $x_{i}(t)$ is the expression of gene $i$, and $P_{i}(t)$ is the
production rate at time $t$. Noise is added through the independent Gaussian
white noise processes $\alpha$ and $\beta$. The production rate is modelled as
a sum of contributions from each parent:

$P_{i}=\sum\limits_{j\in\text{pa}(i)}p_{ij}+b_{i}$ (14)

where $p_{ij}$ is the regulatory effect of parent $j$ on gene $i$ and $b_{i}$
is the basal production rate of gene $i$. $p_{ij}$ is determined by a non-
linear function of the expression of the parent according to:

$p_{ij}=K_{ij}\frac{x_{j}^{n_{ij}}}{h_{ij}^{n_{ij}}+x_{j}^{n_{ij}}}$ (15)

if parent $j$ is an activator and

$p_{ij}=K_{ij}\left(1-\frac{x_{j}^{n_{ij}}}{h_{ij}^{n_{ij}}+x_{j}^{n_{ij}}}\right)$
(16)

if parent $j$ is a repressor. Here, $K_{ij}$ denotes the maximum contribution
of regulator $j$ to target $i$, $n_{ij}$ is the Hill coefficient that
modulates non-linearity, and $h_{ij}$ is the half-maximal effect.

The expression of gene $i$ is updated according to:

$(x_{i})_{t+\Delta t}=(x_{i})_{t}+(P_{i}(t)-\lambda_{i}x_{i}(t))\Delta
t+q_{i}\sqrt{P_{i}(t)}\Delta W_{\alpha}+q_{i}\sqrt{\lambda_{i}x_{i}(t)}\Delta
W_{\beta}$ (17)

where $\Delta W_{\alpha}\sim\sqrt{\Delta t}N(0,1),\Delta
W_{\beta}\sim\sqrt{\Delta t}N(0,1)$. See [12] for full details.

Here, we modify SERGIO to perform hard interventions by setting $x_{i}=a$ for
the $i$-th node, as done by other applications of SERGIO [25].

###  Appendix C Negative Binomial Likelihood

As conducted by [17], the Negative-Binomial likelihood used to fit the
observational mean $\mu_{j}^{0}$ and inverse dispersion $\theta_{j}$ for node
$j$ is given by,

$f_{\text{NB}}(x_{j};\mu^{0}_{j},\theta_{j})=\frac{\Gamma(x_{j}+\theta_{j})}{\Gamma(x_{j}+1)\Gamma(\theta_{j})}\left(\frac{\theta_{j}}{\theta_{j}+x_{j}}\right)^{\theta_{j}}\left(\frac{\mu_{j}}{\mu_{j}+\theta_{j}}\right)^{x_{j}}.$
(18)

The log-likelihood $\log{f_{\text{NB}}(x_{j};\mu_{j}^{0},\theta_{j})}$ is
optimised using the L-BFGS-B implementation in SciPy [19].

###  Appendix D Extended Empirical Results

In this section, we present initial benchmarking results of commonly used CSL
methods using CausalRegNet and a description of the computational resources
needed to perform these experiments in Section D.1. We start by describing the
parameter settings used for simulations involving CausalRegnet and SERGIO in
Section D.2. Section D.3 then shows full simulation results. The methodology
used for selecting the 100 cancer-related genes for CausalRegNet is given in
Section D.4. A comparison of the distribution of correlations in the
CausalRegNet synthetic dataset and Replogle as well as for the causal and non-
causal relationships simulated by CausalRegNet are shown in Section D.5. The
distribution of Wasserstein distances between synthetic and real gene
distributions for the true and random gene labels is given in Section D.6.

####  D.1 Computational Resources

CPU hours were used to run all simulations. The most intense experiments
involved generating data using SERGIO for the scalability results in Section
4.2. This involved over 300 compute hours with over 64GB of active memory. The
remaining experiments in Section 4.3 and Section 4.4 were run locally with
basic compute requirements. The experiments involving CSL algorithms also
exceeded over 150 compute hours to run the algorithms on graphs of size 50 and
100 with 16GB of active memory.

####  D.2 Simulation Parameter Settings

For CausalRegNet, individual node means and inverse dispersions are obtained
from fitting a negative binomial distribution to randomly selected genes in
Section D.4. The minimum effect, maximum effect and regulatory adjustment
effect for each node are $\\{\alpha_{j}=2,\beta_{j}=0.1,a_{j}=10\\}$,
respectively, and edge weights are sampled from the range
$w_{ij}\in(-2,-0.5)\cup(0.5,2)$. For SERGIO, we follow the default parameters
specified by [12], which are given in Table 2. The MR is the master regulator
of each causal graph, that is, the node with no parents in each graph. DAGs
are randomly generated according to Table 3 and identical DAGs are used for
both CausalRegNet and SERGIO.

Table 2: Parameters for synthetic data generation with SERGIO.  Interaction Strength | $k_{ij}\sim\text{Unif}(1,5)$
---|---
MR Basal Production | $b_{i}\sim\text{Unif}(5,15)$
Noise Amplitude | $q=1$
Decay | $\lambda=0.8$
Hill Coefficient (degree of non-linearity) | $n_{ij}=2$
Table 3: Parameters for DAG generation. Repetitions | 10
---|---
Graphs | ER-$m$
Nodes | $d\in\\{10,20\\}$
Edges | $m\in\\{2\cdot d,4\cdot d\\}$
Samples Per Setting | $n=1000$

####  D.3 Benchmark Study

Here, we conduct a simulation study to examine the performance of various CSL
algorithms using CausalRegNet as our synthetic data generator.

Our simulation study includes a diverse set of CSL methods. The PC algorithm,
introduced by [26], is a constraint-based method, utilising conditional
independence testing to identify the causal strucuture with observational
data. [2] developed a continuous optimisation scheme on the space of causal
graphs through NO-TEARS. Finally, we include the GSP and IGSP methods based on
traversing the space of sparsest permutations, which were developed by [27]
and [28], respectively. IGSP is an interventional method and thus will be
provided both observational and interventional data. Each method is
implemented according to the settings shown in Table 4.

#####  D.3.1 Setting

We use a normalised Structural Hamming Distance (SHD), which we define as
$\text{SHD}/(m\cdot d)$, false omission rate (FOR) and false discovery rate
(FDR), to evaluate the performance of each method on full graph retrieval.
Both the Erdős-Rényi (ER) and Scale-Free (SF) graph strucures methods are
used. The full parameters for generating DAGs that are used to generate data
with CausalRegNet are shown in 3. For observational methods, data coming from
only the observational setting, denoted by $\mathcal{I}=\emptyset$, is
provided.

Table 4: Parameter settings used for CSL algorithms for the benchmark. PC | `indep_test = kci`
---|---
`stable` = `True`
`alpha = 0.05`
NO-TEARS | `class = mlp`
`max_iter = 20`
GSP |  `ci_test` = `kci`
`alpha = 1e-3`
`depth = 5`
IGSP |  `ci_test` = `kci`
`inv_test = kci`
`alpha = 1e-3`
`depth = 5`

#####  D.3.2 Results

Figure 6(a) shows that the normalised SHD for all methods increases with graph
size, suggesting that the size of the graph significantly increases the
difficulty of CSL. Furthermore, this trend corresponds to decreasing
specificity of each of the CSL algorithms as shown by increasing FDR with
graph size in Figure 6(c). Figure 6(b) shows that as the graph size grows, all
methods are less likely to omit edges incorrectly, supporting our hypothesis
of falling specificity. NO-TEARS consistently outperforms the other methods.

This benchmark highlights that despite the theoretical foundations of methods
in CSL, there exist settings where provably correct assumptions, such as
causal sufficiency, and causal faithfulness are not enough to ensure good
performance. Thus, the need for testing the sensitivity of methods to
violations of modelling assumptions in a context-specific manner is crucial.

(a) SHD performance for various causal graph settings.

(b) FOR performance for various graph settings.

(c) FDR performance for various graph settings.

Figure 6: Results from benchmark with CausalRegNet.

####  D.4 Selecting Cancer Genes

We utilised Entrez Direct tool [29] and Pangaea [30] to parse abstracts from
PubMed articles using the search term "cancer". From these abstracts, we
identified the 100 most commonly referenced cancer-related genes, which were
both perturbed and measured in the Replogle dataset.

####  D.5 Correlation Patterns

Furthermore, Figure 7(a) shows the correlation between 100 genes simulated by
CausalRegNet and their real-world correlation in [8]. The overall
distributions of correlation have similar shape where they are relatively
centered about some mean correlation value. The real data shows a positive
skew, which suggests that there is more upregulation than downregulation
occurring in this particular set of genes than expected. For this simulation,
CausalRegNet assumes a balanced distribution of upregulators ($w_{ij}>0$) and
downregulators ($w_{ij}<0$). Some tweaking of this balance may be needed to
achieve a more comparable distribution of correlations. This skew may also be
due to all genes studied here being cancer-related. Therefore, they may be
more correlated in expression than would be expected with a random subset of
genes taken from [8]. Further work would be needed to verify this.

We also look at the difference in correlation between nodes that have a direct
causal relationship and those that do not. Figure 7(b) shows that there is a
clear relationship between causal links and correlation.

(a) Distribution of correlations between genes in data simulated by
CausalRegNet and real data from [8].

(b) Distribution of correlations between nodes that have a direct causal
correlation and those that do not.

Figure 7: Correlation patterns in data generated by CausalRegNet.

####  D.6 Wasserstein Distance between Real and Synthetic Data

The results shown in Section 4.4 show that simulating data with fitted
negative binomial parameters in CausalRegNet leads to better than random
marginal distribution fit when compared to real data. Although this is true
for the majority of nodes, there can still be outliers when simulating large
graphs as shown in Figure 8. These outliers occur at a constant rate thus for
smaller graphs, a smaller number of poorly fitting nodes will appear.

Figure 8: Distribution of Wasserstein distances


11institutetext: School of Computing Science, Newcastle University

# Twenty-Five Shades of Greycite: Semantics for referencing and preservation

Phillip Lord  Lindsay Marshall

###### Abstract

Semantic publishing can enable richer documents with clearer, computationally
interpretable properties. For this vision to become reality, however, authors
must benefit from this process, so that they are incentivised to add these
semantics. Moreover, the publication process that generates final content must
allow and enable this semantic content. Here we focus on author-led or “grey”
literature, which uses a convenient and simple publication pipeline. We
describe how we have used metadata in articles to enable richer referencing of
these articles and how we have customised the addition of these semantics to
articles. Finally, we describe how we use the same semantics to aid in digital
preservation and non-repudiability of research articles.

##  1 Introduction

The academic publishing industry is changing rapidly, partly as a result of
external changes such as the move to open access, and partly as a final
recognition in the importance of the web. With change comes the opportunity to
add more semantics to publications [1, 2, 3], to increase the computational
component of papers, enabling publication to take its place in the linked data
environment [4].

While within academia, third party publication — where knowledge is given to a
third party to manage the publication process — is commonplace, outside in
many technical disciplines we see direct publication, where the author
publishes work that readers can then directly access. This form of publication
is often called “grey literature” publication — a somewhat derogatory term —
however, it has some significant advantages. It is rapid and places the author
in control, allowing them to innovate in terms of presentation and content[5].
It operates without editorial control from third-party publishing which may
help to overcome the publication bias found in many areas of scientific
publishing. We have previously used a form of grey literature publishing to
publish ontology tutorial material[6, 7]; this has resulted in the release of
useful material which would otherwise probably not have been created, as many
academics regard book chapters as having little purpose[8].

From the perspective of semantic publishing, it has an additional advantage;
the process is often very simple, without additional human intervention
between the author and the final published form. This simplicity means that
semantics added by the author can pass through to the published version with
relative ease. In the process, it is also possible that semantics added by the
author can aid in the authoring process, which we consider of key
importance[9].

However, grey literature publishing lacks some of the formality of third-party
academic publishing; for instance, several organisations provide centralised
collection of bibliographic metadata; we have used this metadata, for
instance, to enable accurate citation of academic literature through the use
of primary identifiers. The lack of a centralised authority for author
published literature, however, prevents this technique from being used for
general URIs. This presents us with a simple research question: are there
enough semantics on the extant web to provide clear bibliographic metadata for
different web pages?

In this paper, we describe two new systems: greycite and kblog-metadata. The
former, addresses the problem of bibliographic metadata, without resorting to
a single central authority, extracting this metadata directly from URI
endpoints. The latter provides more specialised support for generating
appropriate metadata. We describe how these systems support our three steps
doctrine [9], which suggests that semantic metadata must be of value to all
participants in the publishing process including the authors. We also describe
how these systems can impact on another major problem with author-led
publishing: that of archiving and “link-rot”.

##  2 References

Referencing is ubiquitous within scientific and academic literature, to the
extent that it can be considered to be a defining feature. Academics reference
previous work both as a utility to the reader, and as a mechanism for
establishing provenance. However, reference insertion and formatting is
complex to the point of humour[10]; with nearly 3000 citation formats in
common use[11], reversing the process is even harder.

We have previously described the _kcite_ tool which enables automatic
generation of reference lists from a primary identifiers[9]: as described
previously, it is often possible to hide these from the user behind tooling,
so that they do not need to insert primary identifiers by hand[9]. This form
of referencing also has advantages for human and machine consumpution of the
data; the primary identifier, which is also accessible to downstream analysis;
moreover, because the reference is generated as a result of this identifier,
when the author checks the reference, they are also effectively checking the
identifier, which conventionally, the author must check manually at extra cost
to their time.As an ad hoc measure, user feedback from our tool has now
identified a number of primary identifiers (DOIs) with inaccurate metadata,
and one systematic error in the presentation of these identifiers affecting
many institutional repositories[12].

This, however, requires a source of metadata: currently, kcite supports (most)
DOIs, arXiv and PubMed IDs directly, all of which allow metadata harvesting.
Following the development of kcite, our request resulted in both CrossRef and
DataCite – the two most significant DOI registration agencies for academia –
providing metadata in the form that kcite consumes (Citeproc JSON). For
general URIs, unfortunately, there is no centralised authority which can
provide this metadata.

###  2.1 Technical Glossary

Here, we provide a short technical glossary of the tools described, also shown
in Figure 1, as an aid to understanding.

kcite:

A wordpress plugin that generates a reference list for an article from primary
identifiers. Uses a variety of services, including greycite, to resolve
identifiers to bibliographic metadata.

kblog-metadata:

A wordpress plugin that provides flexible presentation of bibliographic
metadata, both computationally and visibly through on-screen widgets.

greycite:

A server which returns bibliographic metadata for any URI, extracted from that
URI for the article resolved by that URI.

Citeproc JSON:

A bibliographic format defined by the Citeproc-js tool.

BibTeX:

a format defined by the BibTeX tool.

##  3 The Greycite System

We initially considered the possibility that kcite could harvest its own
metadata directly. It would have been possible, for instance, for a kcite
installation on one site to return metadata to another, through a REST call,
or as embedded metadata. However, this would have required users to know in
advance which URIs were so enabled, and would have worked with few websites.

To avoid this limitation, we wished to use extant semantics already on the
web; the complexity of this task argued against integration with kcite which
is an end-user tool. Additionally, as a server greycite would usable to more
than one client; in fact, this has proven to be the case, with a third-party
tool, knitcitations which supports dynamic citations in a literate programming
environment for R[13].

Greycite provides bibliographic metadata in a variety of formats on request
about an arbitrary URI; an architectural overview is shown in Figure 1. It
uses a simple REST API to do this, and returns either Citeproc-JS JSON
(directly used by kcite)[14], BibTeX (used by knitcitations, and the kblog-
metadata tool described here). We have additional support for other formats,
including RDF (encoding Dublin Core), RIS, and Wikipedia “cite” markup. It
store the results of metadata extractions, initially for reasons of
efficiency, although this is also valuable for ephemeral sources of
metadata(see Section 4).

Greycite extracts a number of sources of metadata, and uses a scoring scheme
and a set of heuristics to choose between them; we describe these next.

Figure 1: Client server interaction between Greycite and clients

##  4 On how the Web describes itself

To enable referencing, we need five key pieces of bibliographic metadata,
namely:

  * •

Author(s) (A)

  * •

Title (T)

  * •

Date of publication (D)

  * •

“Container” – equivalent to journal, conference or website. (C)

  * •

Canonical Identifier (I)

These are the minimal pieces of metadata used by most referencing styles, and
following standard publication practices. We have now investigated many
sources of web-delivered metadata. These have been discovered in a number of
ways: some were designed for this purpose. Others, where discovered by
inspection of academic websites; some were discovered entirely by chance
(where an authors name was visible on a web page, but not extractable, we
search for all instances of that name, looking for structure). We prioritised
“interesting” websites, for our definition of interesting.

A complete list of all the mechanisms greycite uses for metadata extraction is
shown in Table 1. By itself HTML provides very few of these five pieces of
metadata; only the title is extractable; even here, for most browsers, the
title is displayed publicly, in the browser title bar. As a result, many sites
include the name of the site, often “—” delimited in the title of each page,
which makes this a relatively messy form of data.

We also investigated the use of CoINS metadata; this standard is used by a
number of academic websites, and can be consumed by a some bibliographic
tools111http://ocoins.info. It is an imperfect tool. The standard is rather
confusing to read, the main website describes it as using a NISO 1.0 Context
Object, the link to the specification for which is broken. Different
implementations tend to produce different variations of the same metadata.
More over, CoINS metadata does not necessarily describe the article being
posted; for example, http://researchblogging.org uses CoINS to describe a
secondary article being reviewed. It has a significant advantage, however,
over most metadata specifications which is that it is embeddable in the _body_
of a web page; for hosted websites, authors often do not control the headers
and cannot add elements to it.

The guidelines for inclusion to Google Scholar are somewhat clearer, and
easier to implement, although even here there are common causes for confusion
(citation_author vs citation_authors). This form of metadata is relatively
common on many journal websites, but is not, in our experience, wide-spread
outside academia. More common, is Open Graph Protocol, or OGP222http://ogp.me;
this is a form of RDFa developed as part of the Facebook platform. It is found
on a large number of websites including many common blog platforms, as well as
various news outlets, such as BBC News, which are otherwise hard to cite. The
author list is often not represented in OGP333Including on the OGP website!;
while OGP has the ability to do this, authorial metadata needs to be gathered
from a secondary URI, linked from the main content; this is more complex to
implement, which may explain why it is commonly missing.

Source | Type |  Notes
---|---|---
Atom | TDCAI |  Inferences where article is not present
CoINS | TCDA |  Blocked where identifier does not match location
CEUR-WS | TCDA |  Uses span tags in index files
Dublin Core | TCDA |  Both dc: and dc. recognised
Eprints | TCDA |
EXIF | TDA |  In Progress
FOAF | N/A |  In Progress
GIF | N/A |  In Progress
Google Scholar | TCDA |  Both citation_author and citation_authors. Bepress prefix with bepress
HTML | T |  The “title” tag
Link | N/A |  In Progress
Meta | TCDA |  Common uses recognised
OGP | TCDAI |  Some syntactic variants
OpenLibrary | TCDI |  Preliminary
ORCID | TCA |  Screen Scraping
PDF | TA |  Often fails!
Prism | CD |
RSS | TCDAI |  See Atom
Schema | TD |  In Progress
Scholarly HTML | N/A |  In Progress. Never seen in the wild
ScienceDirect | TCD |  Screen Scraping
Twitter | TCAI |  “Author” is normally a hashtag
URI | D |  Heuristic based on link structure
W3C | TCDAI |  Screen scraping specific for W3C specifications
WorldCat | TDA |  Screen scraping
Table 1: Twenty-Five Sources of Metadata: type indicates the metadata
extractable (Type, Date,Container,Author,Identifier).In progress indicates
that we believe more metadata is present. Screen Scraping means heuristics
based on HTML structure.

Another source of authorial metadata is RSS/Atom feeds. Many common platforms
include a dc:creator tag and this is often the only easily extractable form of
metadata. We do find that generic (admin, blog admin) or personal but informal
(Phil, phillord) user names are fairly common; this is the default behaviour
for many content management systems, and appears to be a conscious choice for
many multi-user sites. Greycite filters some of the more common ones and does
not consider them as valid metadata. We also provide heuristics where articles
are missing; for instance, if all articles in an RSS feed have the same author
and container title, we infer this information for missing articles.

Another commonly missing piece of metadata is date; while it can be found in
RSS/Atom feeds, these are not always present and are _ephemeral_. In contrast
to author or container information, publication date cannot be infered where
articles are missing from metadata on other articles. We apply a heuristic
here in acknowledgement of the fact that many blogs use a date format for
their URI permalinks. In fact of the URIs in greycite, we can mine date
metadata from some 33% of them; while this is not a representative sample, it
does show that heuristics can be surprisingly effective.

Unfortunately, many scientific papers are published in PDF; while we do
attempt to extract metadata from these, greycite is currently not very
effective, so most PDFs appear to contain no extractable metadata; we are
investigating more PDF parsers to attempt to address this problem. In some
cases, we have provided heuristics which work around this difficulty: greycite
will provide metadata for PDFs hosted by CEUR-WS; however, we achieve this by
mining metadata from the HTML files which link to the PDF.

A significant number of websites do not provide any specific metadata that we
were able to discern; interesting and surprising cases include most of the W3C
standards, websites for both the International and Extended Semantic Web
Conferences, and the ORCID webpages. We have a significant number of special
purpose extraction plugins; for instance, from a desire to reference W3C
specifications, we have created a single site plugin which uses a highly ad
hoc screen-scraping technique. Taken together, of the 4000 URI that have been
submitted, Greycite can extract the main four pieces of metadata (TCDA) from
62% of URIs.

##  5 On how the web could describe itself

While Greycite can extract metadata from many different sources, it does
require some support from the content. Unfortunately, for many content
management systems whether this metadata is available or not is dependent on
the local setup; for instance, with WordPress, the presence or absence of many
sources of metadata is theme dependent; the exception to this is data from the
RSS/Atom feeds although even here, the feeds themselves can be disabled at the
theme level444This would generally be considered to be a broken theme, or
through author choice555This would generally be considered to be a broken
author.

We have therefore created a plugin for WordPress to address this need; while
the solution is, of course, specific to WordPress, the use cases that we
address are considerably more general. This plugin, _kblog-metadata_ ,
currently adds metadata in three formats: Google Scholar, OGP and CoINS. The
latter is used by and has been tested with Zotero and similar bibliographic
tools, which is the main reason for its inclusion. Facebook provide an
explicit tool for testing OGP, while Google Scholar do not. By default, kblog-
metadata requires no configuration and uses knowledge directly from the
WordPress container, which provides suitable values for the five pieces of
metadata we require (see Section 4).

While the author can check that their metadata is appearing correctly, through
the use of greycite, this requires them to use a secondary website.
Alternatively, they can link to their article using kcite, which will then
generate a reference on the basis of the metadata; however, this requires
creating new content, to check old. Following our three steps doctrine, we
wished to make the metadata more useful for the authors. We have, therefore,
added “Widget” support, which displays citation information for each page (or
a website as a whole) using the same metadata resolution techniques; this
display both eases the task of checking the metadata, as well as incentivising
the author to do so. The widget also provides a BibTeX download of the
citation. As well as being useful for authors and readers, this has an
additional utility: the BibTeX actually comes from greycite, on the basis of
its metadata extraction. When anything (including robots) access this BibTeX,
Greycite is invoked, and hence becomes aware of the new article.

Although for simple use, the default WordPress data suffices, there are
several uses cases where it does not. Therefore, kblog-metadata provides
authors with the ability to set the metadata independently on an individual
post basis. This fulfils a number of use cases. The most common of these is
for multiple-author posts; WordPress multiple author support is built around
editing rights, rather than authorship. Hence all authors must have WordPress
logins which they otherwise may neither want or need. Kblog-metadata allows
setting authorship lists independently of login rights. Secondly, authors may
also wish to provide an alternative container title. Combined, these two
facilities enable WordPress to operate as an “preprints” server. For example,
http://www.russet.org.uk/blog/2054 resolves to the full text of our paper[9],
which uses both facilities so that the citation appears with three authors,
and “Sepublica 2012” as the container title.

Since, kblog-metadata was released, WordPress also supports “Guest authors”
through the co-authors-plus plugin – which likewise dissociates login rights
from authorship; this provides a much nicer graphical environment for defining
co-authors than kblog-metadata, but comes with an overhead that authors must
be created individually. Kblog-metadata will use metadata from this plugin if
it is installed.

Finally, we have added support for the use of shortcodes to define authorship.
This is very useful when content is being generated outside of the WordPress
environment; for example, on http://bio-ontologies.knowledgeblog.org, most of
the content is generated from Word documents. During publication, we markup
the author names with shortcodes — [author]Phillip Lord[/author]; this markup
passes unmolested through Word’s HTML conversion and is then interpreted by
WordPress. This prevents cut-and-paste errors that would occur if authors had
to be added manually — a significant issue for science where most articles
have many authors. This website also modifies the container title to
distinguish between different years.

##  6 Identifying by Proxy

One significant issue with kcite as a referencing engine is the requirement
for a primary identifier for every item666kcite does allow addition of all
citation metadata within an inline shortcode, although this is intended as a
fallback. Most scientific literature, and any article posted on the web is
likely to have an identifier that kcite can use. However, this causes problems
for two specific types of resource. First, many smaller conferences and
workshops do not publish their literature in a web capable form; in many cases
papers on the web are available as PDF or Postscript only. And even when web
hosted, sites may not add bibliographic metadata. Kblog-metadata provides a
partial solution to these problems: authors can host their articles, and alter
the metadata accordingly as described in Section 5. However, this fails for
work by other authors, whose work cannot be posted without permission. A
similar problem exists for books; while these generally do have a standard
identifier (ISBN) we have not been able to find a publicly available mechanism
to automate the transformation from ISCN to structured bibliographic metadata.

Greycite provides a mechanism to address this difficulty. There are a number
of catalogues available for both scientific literature and books; these often
have a primary URI which can be used as a reference identifier. Greycite
currently supports several sites of this form: WorldCat (http://worldcat.org)
provides URIs for books (as well as other forms of media such as CDs and
DVDs), Mendeley which references journal articles and OpenLibrary
(http://www.openlibrary.org) which also provides URIs for books. In these
cases, references will appear correctly when used in Kcite, showing the source
of metadata which could, in principle, be used to track the original resource.

##  7 Metadata for Preservation

One recurrent issue with author-led publishing is the difficulties associated
with digital preservation; custom and practice means that it is considerably
harder for author-led publications to ensure that work is preserved than
third-party publications; systems such as CLOCKKS or LOCKKS are often just not
accessible to smaller-scale author-led publication.

To address this need, we have integrated greycite with public archiving
efforts, such as the Internet Archive, the UK Web Archive and WebCite. As well
as scanning URIs for metadata, greycite periodically checks these archive
sites, to see if they are available as archives. We use this metadata in a
number of ways.

First, archive sites are available directly from Greycite through a REST API.
Kblog-metadata provides an “archived” widget where it publicly displays this
information; this provides a third-party stamp that the article has actually
been available from the time stated, as well as an ad hoc enforcement of non-
repudiability. If an author changes their own content, the differences with
the archived sites will be clear.

If a site disappears, then these links to the archives will also disappear.
Kblog-metadata also allows readers to download BibTeX files for any (or all)
articles; this metadata comes directly from Greycite and includes links to all
known public archives. Anyone citing an article using this file will therefore
have a reference to archival versions.

Of the services we currently check for archival versions, currently only
WebCite offers on-demand archiving777WebCite is asking for funding on the web,
which is an unfortunate sign. Greycite currently submits any archive with four
(TCDA) piece of metadata to WebCite for archiving. Additionally, greycite
itself stores historical metadata for indeterminate amounts of time, and
therefore constitutes a metadata archive.

##  8 Tracking movement around the web

In addition to the four pieces of bibliographic metadata, greycite collects
one other key piece of knowledge; a canonical URI. Currently, this knowledge
is not represented in many of the formats we harvest. While, CoINS does
provide a field which can be used for this purpose, in practice it is not that
useful: CoINS is used to embed bibliographic metadata into the web, but the
CoINS may not relate to the article in which it is embedded. Open Graph
Protocol data also returns an explicit identifier; in this case, this is about
the article in question. This means Greycite can store a canonical URI for a
particular article, regardless of the URI used to access the article. Again,
and perhaps unexpectedly, RSS/Atom feeds are extremely useful; these carry a
link and explicitly state whether it is a permalink (i.e. canonical) or not.

The presence of a canonical URI makes it possible to track content as it
moves. For instance, it is relatively common for blogs to change their
permalink structure; with WordPress, for instance, existing links are
maintained through the use of a 301 Redirect response. Greycite could
recognise this situation and use the Redirect location as the canonical link;
unfortunately HTTP redirects are used for many different purposes, including
load balancing. Instead, greycite recognises that the URI used to fetch a
request and the stated canonical URI are different and records this fact. For
example, greycite records that http://www.russet.org.uk/blog/2012/02/kcite-
spreads-its-wings/ changed from being canonical to not sometime between April
2012 and Jan 2013 (actually this happened in June 2012).

Currently, greycite returns the canonical URI with requests for both BibTeX or
JSON data; authors of referring documents will therefore will have an recent
link. Although, currently not implemented, we plan to add more explicit
support for this to our kcite client, so that it will display canonical URIs;
again, this supports digital preservation. Articles which refer to URIs which
have ceased to be canonical, would be able to display both the URI to which
the author original referred and the correct canonical reference.

The ability to track articles as they move also opens up a second possibility.
Currently, one main stated advantage of systems like DOIs is the ability to
change the location of a record without necessitating a change in identifier.
A similar system is also available in the form of PURLs (persistent
URLs)888http://purl.oclc.org/. Greycite allocates PURLs for all URIs for which
it can extract all the required piece of metadata. Currently, these redirect
to the last known canonical URI for a given URI; in effect, this means that
PURLs will track URIs for any website that maintains its redirects and
metadata for sufficient time for Greycite to discover this.

##  9 Discussion

As we have previously stated[9], our belief is that semantic metadata, if it
is to be useful at all, must be useful to all the key players in the
publication process; critically, this includes the author. The tools that we
have described here obey this doctrine; we seek to aid and reward the authors
who use good metadata.

Kcite already follows this principle: if links are inaccurate, then the
reference will not format correctly (or at all). As well as errors made during
authoring, we (PWL) have found non-functioning DOIs, as well as one systematic
error in DOI presentation which has resulted in a change to CrossRef display
guidelines[12]. In addition, correct formatting of the references depends on
the metadata being correct. Again, here, we have found DOIs with inaccurate
metadata. With the addition of greycite, this functionality has been extended
to any URI. Authors are very likely to cite themselves. If they do so, they
are now dependent on their own metadata; if the metadata is wrong, then
references will be. This provides an incentive for authors to correct metadata
for their own purposes, simultaneously making everyone’s life
better999slightly. As well as correcting our own websites, use of greycite has
discovered inaccurate metadata in commercial publishing websites.

Greycite is currently unique so suffers from some of the limitations of
centralisation; however, effectively, it is just a cache. The metadata that it
provides is sourced from the distributed resources that are referenced; it can
support multiple installations trivially. Except in the case of ephemeral
metadata, none of these would be privileged. The current implementation of
greycite also provides an initial answer to our question, is there enough
bibliographic metadata on the web to enable citation: a qualified yes. Through
the use of existing metadata schemes and some heuristics, we can discover this
metadata for many websites. An early analysis suggests that greycite can
provide the four key pieces of metadata for around 1% of the web, which
consitutes 100s of millions of URIs; the percentage for “interesting” websites
is much higher, at over 60%. We currently also lack any statistical analysis
on how correct this metadata is; by inspection, the level of correctness
within the $\sim$4000 URIs submitted from 254 independent IP addresses is
high, but this result is biased as we have corrected errors iteratively. For
random URIs, we lack a gold standard, and most are not in English making
inspection hard.

While using metadata to generate references is useful, it is one-step removed.
The author is not supported in discovering that their metadata is inaccurate
until sometime after it has been published. Kblog-metadata now improves on
this process and makes it more immediate; by visualising the metadata on
publication, authors can check that it is correct. Likewise, the same metadata
is used to generate a BibTeX file which they can use. As an open source tool,
it is hard to know how many installations kblog-metadata currently has,
although download statistics would suggest 30 or 40, including one journal.

Set against this desire to improve the quality of metadata on the web,
greycite has taken a pragmatic approach to the metadata standards it uses. It
currently supports many of the different ways of marking up bibliographic
metadata. More over, it uses many heuristics, to cope with metadata which is
unclean or just broken. This works against the notion of encouraging authors
to improve their metadata; however, increasing the utility of the API makes
this a compromise well worth making.

We are also addressing the issue of digital preservation; we achieve this in
two ways. First, we leverage existing web archives, deep linking through to
them where content has already been archived. To achieve this in a simple
manner requires no semantics at all, beyond the URI for a given resource.
However, resources may be present at more than one URI, or may change their
canonical URI over time. Greycite is now making preliminary use of this
metadata to track articles as they move; the current location can be retrieved
by a client, or alternatively greycite provides PURLs which will work with any
client.

As with our previous work, the level of semantics provided or used by these
publication tools is not high; however, by using existing metadata standards,
greycite can provide metadata for 100s of millions of URIs including many from
websites which are unlikely to care about academic referencing. We have
focused on adding value for authors, both when referencing or displaying
citations on an article. By adding value for the authors, we help to ensure
that they will add value to the metadata. While this approach adds very small
amounts of metadata for an individual article, the aggregate total of metadata
over all articles is, potentially, vast.

## References

  * [1] Shadbolt, N., Hall, W., Berners-Lee, T.:  The semantic web revisited.  Intelligent Systems, IEEE 21(3) (2006) 96–101 http://eprints.soton.ac.uk/262614/.
  * [2] Shotton, D.:  Semantic publishing: the coming revolution in scientific journal publishing.  Learned Publishing 22(2) (2009) 85–94 http://delos.zoo.ox.ac.uk/pub/2009/publications/Shotton˙Semantic˙publis%hing˙evaluation.pdf.
  * [3] Shotton, D., Portwin, K., Klyne, G., Miles, A.:  Adventures in semantic publishing: exemplar semantic enhancements of a research article.  PLoS computational biology 5(4) (2009) e1000361 http://dx.doi.org/10.1371/journal.pcbi.1000361.
  * [4] Bizer, C., Heath, T., Berners-Lee, T.:  Linked data-the story so far.  International Journal on Semantic Web and Information Systems (IJSWIS) 5(3) (2009) 1–22 http://eprints.soton.ac.uk/271285/.
  * [5] Zhu, Y., Procter, R.:  Use of blogs, Twitter and Facebook by PhD Students for Scholarly Communication: A UK study.  In: China New Media Communication Association Annual Conference. (2012) https://www.escholar.manchester.ac.uk/item/?pid=uk-ac-man-scw:187789.
  * [6] Lord, P.:  Ontogenesis.  http://ontogenesis.knowledgeblog.org/ (2012)
  * [7] Lord, P., Cockell, S., Swan, D.C., Stevens, R.:  Ontogenesis knowledgeblog: Lightweight semantic publishing.  http://www.russet.org.uk/blog/1920 (2011)
  * [8] Bishop, D.:  How to bury your academic writing.  http://deevybee.blogspot.com/2012/08/how-to-bury-your-academic-writing.%html (2012)
  * [9] Lord, P., Cockell, S., Stevens, R.:  Three steps to heaven.  http://www.russet.org.uk/blog/2054 (2012)
  * [10] Psyphago, .:  Scientists receive 12.6 million dollar grant to format references correctly.  http://collectivelyunconscious.wordpress.com/2013/02/27/scientists-rece%ive-12-6-million-dollar-grant-to-format-references-correctly/ (2013)
  * [11] Zelle, R.:  Styles.  http://citationstyles.org/styles/
  * [12] Lord, P.:  The evil a space can do.  http://www.russet.org.uk/blog/2340 (2013)
  * [13] Boettiger, C.:  Semantic citations for the notebook and knitr.  http://www.carlboettiger.info/2013/02/22/semantic-citations-for-the-not%ebook-and-knitr.html (2013)
  * [14] Bennet, F.:  The citeproc-js processor (2009) https://bitbucket.org/fbennett/citeproc-js.