\section{Proof of Theorem \ref{theorem:approachEquality}}

\input{critical_proof.tex}

\section{Proof of Theorem \ref{newTheoremSimple}}

\input{score_proof.tex}

\onecolumn
\section{Extended Example of Critical Instance}
Here is the complete listing of the critical instance of schema $S_1$, and rule $r_2$ from our example.

\footnotesize \begin{lstlisting}

$\mathds{C}(S_1,r_2)$ = $\{$
 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$,
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$

 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,

 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,

 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\ensuremath{{:}{\lambda}}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\ensuremath{{:}{\lambda}}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\uri{CO\_Danger}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\uri{WorkerTag}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\uri{TunnelA}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}$
 $\gpattern{\texttt{sosa:observedProperty}}{\texttt{sosa:hasResult}}{\lit{1}}$

 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\ensuremath{{:}{\lambda}}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\uri{CO\_Danger}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\uri{WorkerTag}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\uri{TunnelA}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}$
 $\gpattern{\uri{CO\_Danger}}{\texttt{sosa:hasResult}}{\lit{1}}$

 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\ensuremath{{:}{\lambda}}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\uri{CO\_Danger}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\uri{WorkerTag}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\uri{TunnelA}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}$
 $\gpattern{\uri{WorkerTag}}{\texttt{sosa:hasResult}}{\lit{1}}$

 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\ensuremath{{:}{\lambda}}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\uri{CO\_Danger}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\uri{WorkerTag}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\uri{TunnelA}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}$
 $\gpattern{\texttt{sosa:hasFeatureOfInterest}}{\texttt{sosa:hasResult}}{\lit{1}}$

 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\ensuremath{{:}{\lambda}}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\uri{CO\_Danger}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\uri{WorkerTag}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\uri{TunnelA}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}$
 $\gpattern{\uri{TunnelA}}{\texttt{sosa:hasResult}}{\lit{1}}$

 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\ensuremath{{:}{\lambda}}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\uri{CO\_Danger}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\uri{WorkerTag}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\uri{TunnelA}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}$
 $\gpattern{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}{\lit{1}}$

 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\texttt{sosa:observedProperty}}$
 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\uri{CO\_Danger}}$
 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\uri{WorkerTag}}$
 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasFeatureOfInterest}}$
 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\uri{TunnelA}}$
 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\texttt{sosa:hasResult}}$
 $\gpattern{\ensuremath{{:}{\lambda}}}{\texttt{sosa:hasResult}}{\lit{1}}$
$\}$
\end{lstlisting}

\noindent The critical instance of $\mathds{C}(S_1,r_1)$ would be of similar length, and it would feature the \texttt{URI}   \uri{OffLimitArea} instead of the literal \lit{1}.

\section{Introduction}

Inference rules are a common tool in many areas where they are used, for example, to model access control policies \cite{Beimel2011OWLandSWRLAccessControl} and business rules \cite{Fortineau2014SWRLapplicationBusinessRules}. In this paper we are motivated by their use in Internet of Things (IoT) applications, where rules are often used to capture human decision making in a simple and straightforward way \cite{Perera2014RulesReview}. This is especially true in safety-critical domains, such as in the field of Occupational Health and Safety (OHS). OHS knowledge is codified by experts into policies, which are then translated into rules to monitor and regulate workplaces. For example, OHS regulations set limits on human exposure to certain gases. Monitoring systems can use these rules to determine, from sensor data, whether dangerous gas concentration levels have been reached, and trigger warnings or perform actions such as increasing ventilation. Real-world use cases have been provided by industrial partners, such as a wireless and electronic solutions company, and OHS policies from the International Labour Organisation \cite{ILO6331}.

An important limitation of current inference rule applications in the IoT domain is that they require expert human interventions not only to create rules, but also to manage them. This includes determining when they are applicable and what types of facts we can ultimately infer. In the gas-concentration example above, human experts would be employed to answer questions such as: could the rule that aggregates sensor data be used, maybe in conjunction with other rules, to determine whether an area should be evacuated? Is this rule applicable to the available data sources? And will this rule still be applicable after the underlying data sources change (e.g., in case sensors stop working or are substituted with others)? Knowing which rules can be applied and what type of facts can be inferred on a dataset can have safety critical implications, as OHS policies might depend on the availability of certain pieces of information. It is important to note that by executing rules on a specific dataset we only discover the facts that are currently entailed. To predict which facts could potentially be inferred in future versions of the dataset, we need to reason about its schema.

As IoT scenarios become increasingly complex and dynamic, managing rules in a timely and cost effective way requires improvements in automation. In this paper we present an approach that can answer these questions automatically, by reasoning about an abstraction of the available data sources, called the \emph{triplestore schema}. We define triplestore schemas as abstract signatures of underlying data, similar to database schemas. Triplestore schemas can be defined by experts or, as we will see later, can be derived from the types of sensor available. Such schemas are \emph{dynamic}, changing as the data sources change; e.g., a new sensor is added to the system creating triples with new predicates, and entailing other facts.

We consider \texttt{RDF} \cite{Hayes2014GeneralisedRDF} triplestores and we model triplestore schemas as sets of \texttt{SPARQL} \cite{2008SPARQL} \emph{triple patterns}, which in some formal sense restrict or model the underlying \texttt{RDF} data. We express rules as \texttt{SPARQL} \texttt{construct} queries. This type of rules model \texttt{SPIN} \cite{knublauch2011spin} inference rules, which correspond to Datalog rules \cite{Ceri1989Datalog}, and are also compatible with the monotonic subsets of other rule languages, such as \texttt{SWRL} \cite{horrocks2004swrl}. Given an input triplestore schema and a set of rules, our objective is to decide whether the rules would apply on some \texttt{RDF} dataset modelled by this schema. We do so by computing the ``output'' or ``consequence'' schema of these hypothetical rule applications: this is the schema that models all possible \texttt{RDF} datasets that can be obtained by executing the rules on all datasets modeled by the input schema. It is worth noting that our approach is only concerned with schemas, and it is compatible with relational datasets, as long as their schema and and rules can be expressed using \texttt{RDF} and \texttt{SPARQL} \cite{Calvanese2017Ontop}.

Reasoning at the schema level has been explored previously for databases \cite{marnette2009generalized} and Description Logics \cite{glimm2014abstraction}. In fact, for a different but closely related problem of reasoning on database schemas (called \emph{chase termination}), Marnette \cite{marnette2009generalized} employed a canonical database instance, called the \emph{critical instance}, which is representative of all database instances of the given schema, on which we base one of our solutions.

We propose two approaches to reason about the applicability of inference rules on triplestore schemas. First, we re-use the critical instance for our triplestore schemas and develop an approach based on this representative \texttt{RDF} graph: running the rules on this graph produces evaluation mappings which, after careful manipulation in order to account for peculiarities of \texttt{RDF} literals, help produce our consequence schemas. When constructing the critical instance, as in the original case of relational databases, we need to place all constants appearing in our schema and rules in the constructed instance in many possible ways.
This leads to a blowup in its size and so we turn our attention to finding a much smaller representative \texttt{RDF} graph, that we call the \emph{sandbox} graph, and which we populate with only one ``representative'' element. We then develop a novel query rewriting algorithm that can compute the consequence schema on the sandbox graph. We provide correctness, complexity, and evaluation results and experimentally exhibit the efficiency and scalability of our rewriting-based approach: it surpasses the critical-instance methods by orders of magnitude while scaling to hundreds of rules and schema triples in times ranging from milliseconds to seconds.

\section{Background} \label{background}

We consider triplestores containing a single \texttt{RDF} \emph{graph}, without blank nodes. Such a graph is a set of \emph{triples} $\mathds{U} \times \mathds{U} \times (\mathds{U} \cup \mathds{L})$ where $\mathds{U}$ is the set of all \texttt{URI}s, $\mathds{L}$ the set of all literals and $\mathds{V}$ the set of all variables. We use the term \emph{constants} to refer to both literals and \texttt{URI}s. A \emph{graph pattern} is a set of \emph{triple patterns} defined in: $(\mathds{U} \cup \mathds{V}) \times (\mathds{U} \cup \mathds{V}) \times (\mathds{U} \cup \mathds{L} \cup \mathds{V})$.
Given a pattern $P$, $vars(P)$ and $const(P)$ are the sets of variables and constants in the elements of $P$, respectively. We represent \texttt{URI}s as namespace-prefixed strings of characters, where a namespace prefix is a sequence of zero or more characters followed by a column e.g.\ \uri{a}; literals as strings of characters enclosed in double-quotes, e.g.\ \lit{l}, and variables as strings of characters prefixed by a question-mark, e.g.\ \var{v}. The first, second and third elements of a triple $t$ are called, respectively, \emph{subject}, \emph{predicate} and \emph{object}, and are denoted by $t[x]$, $x \in \tau$ with $\tau$ denoting throughout the paper indexes $\tau = \{1,2,3\}$.

A \emph{variable substitution} is a partial function $\mathds{V} \pfun \mathds{V} \cup \mathds{U} \cup \mathds{L}$.
A \emph{mapping} is a variable substitution defined as $\mathds{V} \pfun \mathds{U} \cup \mathds{L}$. Given a mapping $m$, if $m(?v) = n$, then we say $m$ contains \emph{binding} $?v \to n$. The domain of a mapping $m$ is the set of variables $dom(m)$.
Given a triple or a graph pattern $p$ and a variable substitution $m$ we abuse notation and denote by $m(p)$ the pattern generated by substituting every occurrence of a variable $?v$ in $p$ with $m(?v)$ if $?v\in dom(m)$ (otherwise $?v$ remains unchanged in $m(p)$).

Given a graph pattern $P$ and a graph $G$, the \texttt{SPARQL} evaluation of $P$ over $G$, denoted with $\llbracket P\rrbracket_G$, is a set of mappings as defined in \cite{Perez2006Semantics}.
 A graph pattern \emph{matches} a graph if its evaluation on the graph returns a non-empty set of mappings.
We consider inference rules $A \rightarrow C$, where $A$ and $C$ are graph patterns, and can be expressed as \texttt{SPARQL} \texttt{construct} queries. Note that essentially both $A$ and $C$ in a rule are conjunctive queries~\cite{abiteboulbook}. The \emph{consequent} $C$ of the rule is represented in the \texttt{construct} clause of the query, which is instantiated using the bindings obtained by evaluating the \emph{antecedent} $A$, expressed in the \texttt{where} clause. A single application of a rule $r: A \rightarrow C$ to a dataset $I$, denoted by $r(I)$, is
$I \cup \bigcup_{m \in \llbracket A \rrbracket_I} \left \{ m(C), \; \text{if} \; m(C)\;  \text{is a valid \texttt{RDF} a graph} \right \} $.
Rule notations such as \texttt{SPIN} and \texttt{SWRL} can be represented in this format \cite{Nick2018SWRL2SPIN}. The closure, or saturation, of a dataset $I$ under a set of inference rules $R$, denoted by $clos(I,R)$, is the unique dataset obtained by repeatedly applying all the rules in $R$ until no new statement is inferred, that is,  $clos(I,R) = \bigcup_{i=0}^{i=\infty} I_i$, with  $I_0 = I$, and $I_{i+1} = \bigcup_{r\in R}\{r(I_i)\}$.

\section{Problem Description} \label{sec:probdescr}

To reason about schemas we need a simple language to model, and at the same time restrict, the type of triples that an \texttt{RDF} graph can contain. It should be noted that, despite the similarity in the name, the existing \texttt{RDF} \emph{schema} (\texttt{RDFS}) vocabulary is used to describe ontological properties of the data, but not designed to restrict the type of triples allowed in the dataset.
In this paper we define a \emph{triplestore schema} (or just \emph{schema}) $S$ as a pair $\langle S^G, S^\Delta \rangle$, where $S^G$ is a set of triple patterns, and $S^\Delta$ is a subset of the variables in $S^G$ which we call the \emph{no-literal} set. Intuitively, $S^G$ defines the type of triples allowed in a database, where variables act as wildcards, which can be instantiated with any constant element.

To account for the restrictions imposed by the \texttt{RDF} data model, the no-literal set $S^\Delta$ defines which variables cannot be instantiated with literals, thus $S^\Delta$ must at least include all variables that occur in the subject or predicate position in $S^G$. For example, if $\trip{\var{v1}}{\uri{a}}{\var{v2}} \in S^G_{'}$ and $\var{v2} \not \in S^\Delta_{'}$, then the instances of schema $S_{'}$ can contain any triple that has \uri{a} as a predicate. If $\trip{\uri{b}}{\uri{c}}{\var{v3}} \in S^G_{'}$ and $\var{v3} \in S^\Delta_{'}$, the instances of $S_{'}$ can contain any triple that has \uri{b} as a subject, \uri{c} as a predicate, and a \texttt{URI} as an object. To prevent the occurrence of complex interdependencies between variables, we restrict each variable to occur only once (both across triples, and within each triple) in $S^G$ and in the rule consequents.

A graph $I$ \emph{is an instance of a schema} $S$ if
for every triple $t^I$ in $I$ there exists a triple pattern $t^S$ in $S^G$, and a mapping $m$ such that (1) $m(t^S) = t^I$ and (2) $m$ does not bind any variable in $S^\Delta$ to a literal. In this case we say that $S$ \emph{models} graph $I$ (and that each triple $t^S$  models triple $t^I$). All instances of $S$ are denoted by $\mathds{I}(S)$. We say that two schemas $S$ and $S'$ are semantically equivalent if they model the same set of instances (formally, if $\mathds{I}(S) = \mathds{I}(S')$).
Notice that any subset of an instance of a schema is still an instance of that schema.
A rule $r$ within a set of rules $R$ is \emph{applicable} with respect to a triplestore schema $S$ if there exists a graph $I$ instance of $S$, such that the precondition of $r$ matches $clos(I,R-r)$.

Consider the following example scenario of a mine where sensors produce data modelled according to the Semantic Sensor Network Ontology (SSN) \cite{Taylor2017SSN},
with namespace \texttt{sosa}.
In SNN, sensor measurements are called \emph{observations}. A simple approach to create the schema of a dataset is the following. A dataset that can be populated with sensor measurements of a property \uri{x} (e.g., temperature) can be defined with triple pattern $\gpattern{\var{v1}}{\texttt{sosa:observedProperty}}{\uri{x}}$. Pattern $\gpattern{\var{v1}}{\texttt{sosa:hasResult}}{\var{v2}}$ indicates that the results of these measurements are collected and pattern $\gpattern{\var{v1}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{y}}$ indicates that the results are applicable to a specific entity \uri{y} (e.g., a room).
Similar patterns are presented in \cite{Chaochaisit2016CSVX} in the context of converting \texttt{CSV} data into the SNN format .
In this example, the sensors are deployed only in one tunnel, namely $\uri{TunnelA}$, and schema $S_1$ is:
\scriptsize \begin{lstlisting}
$S^G_1$ = $\{ \gpattern{\var{v1}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$,
      $\gpattern{\var{v2}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
      $\gpattern{\var{v3}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
      $\gpattern{\var{v5}}{\texttt{sosa:hasResult}}{\var{v4}}\}$
$S^\Delta_1$ = $\{\var{v1}, \var{v2}, \var{v3}, \var{v5}\}$
\end{lstlisting} \normalsize

We now consider instance $I_1$ of schema $S_1$. In this instance, the sensors in tunnel A observed both a dangerous gas concentration, represented by the value \lit{1}, and the presence of worker \uri{John}.
\scriptsize \begin{lstlisting}
$I_1$ = $\{ \trip{\uri{o1}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$,
      $\trip{\uri{o1}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
      $\trip{\uri{o1}}{\texttt{sosa:hasResult}}{``1"}$,
      $\trip{\uri{o2}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
      $\trip{\uri{o2}}{\texttt{sosa:hasFeatureOfInterest}}{\uri{TunnelA}}$,
      $\trip{\uri{o2}}{\texttt{sosa:hasResult}}{\uri{John}}\}$
\end{lstlisting} \normalsize

Consider two rules $r_1$ and $r_2$. The first one detects when workers trespass on an ``off-limit'' area, and the second one labels areas with dangerous gas concentrations as ``off-limit''.

\scriptsize \begin{lstlisting}
$r_1$ = $\{ \gpattern{\var{v1}}{\texttt{sosa:observedProperty}}{\uri{WorkerTag}}$,
      $\gpattern{\var{v1}}{\texttt{sosa:hasFeatureOfInterest}}{\var{v2}}$,
      $\gpattern{\var{v1}}{\texttt{sosa:hasResult}}{\var{v3}}$,
      $\gpattern{\var{v2}}{\texttt{rdf:type}}{\uri{OffLimitArea}}\}$
      $\rightarrow$ $\{ \gpattern{\var{v2}}{\texttt{rdf:type}}{\uri{TrespassedArea}}\}$
$r_2$ = $\{ \gpattern{\var{v1}}{\texttt{sosa:observedProperty}}{\uri{CO\_Danger}}$,
      $\gpattern{\var{v1}}{\texttt{sosa:hasFeatureOfInterest}}{\var{v2}}$,
      $\gpattern{\var{v1}}{\texttt{sosa:hasResult}}{``1"}\}$
      $\rightarrow$ $\{ \gpattern{\var{v2}}{\texttt{rdf:type}}{\uri{OffLimitArea}}\}$
\end{lstlisting} \normalsize

Since the precondition of rule $r_2$ matches dataset $I_1$, we can apply the rule and derive a new fact: $\gpattern{\uri{TunnelA}}{\texttt{rdf:type}}{\uri{OffLimitArea}}$. On the instance extended by this new fact, rule $r_1$ is applicable and adds $\gpattern{\uri{TunnelA}}{\texttt{rdf:type}}{\uri{TrespassedArea}}$.

Our approach relies on being able to decide which rules are applicable on a specific triple store schema, e.g., $S_1$, in absence of any particular instance, e.g., $I_1$. Since the precondition of rule $r_2$ matches dataset $I_1$, this rule is directly applicable on schema $S_1$, and we would like to be able to decide this by only looking at $S_1$. Moreover if we can decide this and extend schema $S_1$ with a triple pattern that is the schema of  $\gpattern{\uri{TunnelA}}{\texttt{rdf:type}}{\uri{OffLimitArea}}$ (in this case that schema would be the same triple itself), then we would able to reason with this new schema and decide that rule $r_1$ is also applicable.
In practice, what we would like to do is to compute a schema that captures all consequences of applying our set of rules on any potential instance.

The following definition captures this intuition.
Given a schema $S$ and a set of rules $R$, a schema $S^{\prime}$ is a schema consequence of $S$ with respect to $R$, denoted $con(S,R)$, if $\mathds{I}(S^{\prime}) = \bigcup_{I \in \mathds{I}(S)} \{I' | I' \subseteq clos(I,R)$\}. We can notice that since every subset of an instance of a schema is still an instance of that schema, a dataset can contain the consequence of a rule application without containing a set of triples matching the antecedent. This situation is commonly encountered when some triples are deleted after an inference is made.

Keeping track of the schema consequences allows us to directly see which rules are applicable to instances of a schema without running the rules on the data.
In correspondence to a single rule application $r(I)$, of a rule $r$ on an instance $I$, we define a \emph{basic consequence} of a schema $S$ by a rule $r$, denoted by $r(S)$, as a finite schema $S^{\prime}$ for which $\mathds{I}(S') = \bigcup_{I \in \mathds{I}(S)} \{I' | I' \subseteq r(I)\}$. It is now easy to see that the consequence schema for a set of rules $con(S,R)$ is obtained by repeatedly executing $r(S)$ for all $r \in R$ until no new pattern is inferred. Formally, $con(S,R) = \bigcup_{i=0}^{i=n} S_i$, with  $S_0 = S$, and $S_{i+1} = \bigcup_{r\in R}\{r(S_i)\}$, and $S_n= S_{n-1}$ (modulo variable names).
In the following section we focus on the core of our problem which is computing a single basic schema consequence $r(S)$, and describe two approaches for this, namely Schema Consequence by Critical Instance ($\ensuremath{\texttt{critical}}(S,r)$), and Schema Consequence by Query Rewriting ($\ensuremath{\texttt{score}}(S,r)$).

\section{Computing the Basic Schema Consequence}
\label{sec:approach}

\input{critical.tex}

Having defined how the \ensuremath{\texttt{critical}}\ and \ensuremath{\texttt{score}}\ approaches compute a set of mappings, we now describe the details of the last two phases required to compute a basic schema consequence.

\noindent {\bf Filtering of the mappings} This phase deals with processing the mappings computed by either  \ensuremath{\texttt{critical}}\ or \ensuremath{\texttt{score}}, namely $\llbracket A\rrbracket_{\mathds{C}(S,r)}$ or $\llbracket\mathds{Q}(A)\rrbracket_{\mathds{S}(S)}$. It should be noted that it is not possible to simply apply the resulting mappings on the consequent of the rule, as such mappings might map a variable in the subject or predicate position to a literal, thus generating an invalid triple pattern. Moreover, it is necessary to determine which variables should be included in the no-literal set of the basic schema consequence.
The schema $S^{\prime}$, output of our approaches, is initialised with the same graph and no-literal set as $S$ (i.e.\ $S^{\prime G} = S^G$, $S^{\prime \Delta} = S^\Delta$). We then incrementally extend $S^{\prime}$ on a mapping-by-mapping basis until all the mappings have been considered, at which point, $S^{\prime}$ represents the final output of our approach.

For each mapping $m$ in $\llbracket A\rrbracket_{\mathds{C}(S,r)}$ or $\llbracket\mathds{Q}(A)\rrbracket_{\mathds{S}(S)}$, we do the following. We create a temporary no-literal set $\Delta^m$. This set will be used to keep track of which variables could not be bound to any literals if we evaluated our rule antecedent $A$ on the instances of $S$, or when instantiating the consequence of the rule. We initialise $\Delta^m$ with all the variables of our rule $A\rightarrow C$ that occur in the subject or predicate position in some triple of $A$ or $C$, as we know that they cannot be matched to or instantiated with literals.

Then, we consider the elements that occur in the object position in the triples $t_A$ of $A$. We take all the rewritings $t_q$ of $t_A$ in $\mathds{Q}(A)$ (if using \ensuremath{\texttt{critical}} , it would be enough to consider a single rewriting $t_q$ with $t_q = t_A$). Since the mapping $m$ has been computed over the canonical instance ($\mathds{S}(S)$ or $\mathds{C}(S,r)$ depending on the approach), we know that there exists at least one $t_q$ such that $m(t_q)$ belongs to the canonical instance.
We compute the set of schema triples $t^S$ that model $m(t_q)$, for any of the above $t_q$. Intuitively, these are the schema triples that enable $t_A$, or one of its rewritings, to match the canonical instance with mapping $m$. If $t_A[3]$ is a literal $l$, or a variable mapped to a literal $l$ by $m$, we check if there exists any $t^S$ from the above such that $t^S[3] = l$ or $t^S[3]$ is a variable that allows literals (not in $S^\Delta$). If such triple pattern doesn't exist, then $m(A)$ cannot be an instance of $S$ since it has a literal in an non-allowed positions, and therefore we filter out or \emph{disregard} $m$. If $t_A[3]$ is a variable mapped to $\ensuremath{{:}{\lambda}}$ in $m$, we check whether there exists a $t^S$ such that $t^S[3]$ is a variable that allows literals (not in $S^\Delta$).
If such $t^S$ cannot be found, we add variable $t_A[3]$ to $\Delta^m$. Intuitively, this models the fact that $t_A[3]$ could not have been bound to literal elements under this mapping.
Having considered all the triples $t_A \in A$ we filter out mapping $m$ if it binds any variable in $\Delta^m$ to a literal. If $m$ hasn't been filtered out, we say that rule $r$ is applicable, and we use $m$ to expand $S^{\prime}$.

\begin{figure}[!t]

        \center{\subfloat[]{%
        \includegraphics[clip,width=0.9\columnwidth]{GPPGvsCriticalComparison1Fixb.png}%
        }}\\[-0.1ex]

        \subfloat[]{%

        \includegraphics[clip,width=0.9\columnwidth]{GPPGvsCriticalComparison2Fix.png}%

        }\\[-1.2ex]

        \caption{\label{fig:my-label} Average time to compute 20 schema consequences using \ensuremath{\texttt{score}}\ and \ensuremath{\texttt{critical}}\ as the schema size $\lvert S\rvert$ grows.
        The other parameters are: $\lvert P\rvert=1.5\lvert S\rvert$, $\pi_C=0.1$, $\lvert U\rvert =\lvert L\rvert =\lvert S\rvert$, $\lvert R\rvert =4$, $n_A=2$. Due to large difference in performance, subplots (a) and (b) focus, respectively, on \ensuremath{\texttt{critical}}\ and \ensuremath{\texttt{score}} .}
\end{figure}

\noindent {\bf Schema Expansion.}
For each mapping $m$ that is not filtered out, we compute the substitution $s^m$, which contains all the bindings in $m$ that map a variable to a value other than $\ensuremath{{:}{\lambda}}$, and for every binding $?v \to \ensuremath{{:}{\lambda}}$ in $m$, a variable substitution $?v \to ?v^*$ where $?v^*$ is a fresh new variable. We then add triple patterns $s^m(m(C))$ to $S^{\prime G}$ and then add the variables $s^m(\Delta^m) \cap vars(S^{\prime G})$ to $S^{\prime \Delta}$.

Although the schema consequences produced by $\ensuremath{\texttt{score}}(S,r)$ and $\ensuremath{\texttt{critical}}(S,r)$ might not be identical, they are semantically equivalent. This notion of equivalence is captured by the following theorem.
\input{theorem_equality.tex}

The \ensuremath{\texttt{score}}\ approach (and by extension also \ensuremath{\texttt{critical}} , by Theorem \ref{theorem:approachEquality}) is sound and complete. The following theorem captures this notion by stating the semantic equivalence of $\ensuremath{\texttt{score}}(S,r)$ and $r(S)$.
\input{theorem_completeness_soundness_simple.tex}
For our proofs, we refer the reader to our appendix.

\noindent{\bf{Termination.}}  It is easy to see that our approaches terminate since our rules do not contain existential variables, and do not generate new \texttt{URI}s or literals (but just fresh variable names). After a finite number of iterations, either approach will only generate isomorphic (and thus equivalent) triple patterns.

\noindent{\bf{Complexity.}} Our central problem in this paper, computing the schema consequence for a set of rules, can be seen as a form of datalog evaluation~\cite{abiteboulbook} on our critical or sandbox instance. Datalog has been extensively studied in databases and its corresponding evaluation decision problem (whether a given tuple is in the answer of a datalog program) is known to be EXPTIME-complete in the so called combined complexity~\cite{vardi82}, and PTIME-complete in data complexity ~\cite{vardi82,dantsin2001complexity}. Data complexity in databases refers to the setting in which the query (or datalog program) is considered fixed and only the data is considered an input to the problem.
In our setting, data complexity refers to the expectation that the overall number and complexity of the rules remains small.
It is not hard to see that the corresponding decision problem, stated below, remains PTIME-complete in data complexity. The intuition is that once we construct the critical instance in polynomial time (or alternatively, we grow our set of rules by a polynomial rule rewriting) we have essentially an equivalent problem to datalog evaluation (our rules being the datalog program, and the canonical instance being the database).
\begin{customthm}{3}
Given triple pattern $t_{S^{\prime}}$
and schema $S$ as inputs, the problem of deciding whether $t_{S^{\prime}}$ is in the consequence schema of $S$ for a fixed set of rules $R$ is PTIME-complete.
\end{customthm}

\section{Experimental Evaluation}

\begin{figure}[!t]

         \center{\includegraphics[clip,width=\columnwidth]{Scalability2.png}}

        \caption{\label{figure:scalabilityScore} Average time to compute 20 schema consequences using \ensuremath{\texttt{score}}\ as the number of rules $\lvert R\rvert$ grows, for three configuration of $n_A$.
        The other parameters are set as follows: $\lvert S\rvert =50$, $\lvert P\rvert=60$, $\pi_C =0.1$, $\lvert U\rvert =\lvert L\rvert =\lvert S\rvert$. }
\end{figure}

We developed a Java implementation of the \ensuremath{\texttt{score}}\ and \ensuremath{\texttt{critical}}\ approaches and evaluated them on synthetic datasets to compare their scalability.
We developed a synthetic schema and rule generator that is configurable with 7 parameters: $\pi_C,\lvert P\rvert,\lvert U\rvert,\lvert L\rvert,\lvert S\rvert,\lvert R\rvert, n_A$,
which we now describe. To reflect the fact that triple predicates are typically defined in vocabularies, our generator does not consider variables in the predicate position. Random triple patterns are created as follows. Predicate \texttt{URI}s are randomly selected from a set of \texttt{URI}s $P$. Elements in the subject and object position are instantiated as constants with probability $\pi_C$, or else as new variables. Constants in the subject positions are instantiated with a random \texttt{URI}, and constants in the object position with a random \texttt{URI} with $50\%$ probability, or otherwise with a random literal. Random \texttt{URI}s and literals are selected, respectively, from sets $U$ and $L$ ($U \cap P = \emptyset$). We consider chain rules where the triples in the antecedent join each other to form a list where the object of a triple is the same as the subject of the next. The consequent of each rule is a triple having the subject of the first triple in the antecedent as a subject, and the object of the last triple as object. An example of such rule generated by our experiment is:
$\{ \trip{\var{v0}}{\uri{m1}}{\var{v1}}, \trip{\var{v1}}{\uri{m3}}{\var{v2}} \} \rightarrow \{\trip{\var{v0}}{\uri{m2}}{\var{v2}}\} $
In each run of the experiment we populate a schema $S$ and a set of rules $R$ having $n_A$ triples in the antecedent. To ensure that some rules in each set are applicable, half of the schema is initialized with the antecedents triples of randomly selected rules. The other half is populated with random triple patterns.
We initialize $S^\Delta$ with all the variables in the subject and predicate position in the triples of $S$.  The code used in this experiments is avaliable on GitHub;\footnote{\url{https://github.com/paolo7/ap2}} it uses Apache Jena\footnote{\url{https://jena.apache.org/}} to handle query execution. We run the experiments on a standard Java virtual machine running on an Ubuntu 16.04 computer with 15.5 GB RAM, an Intel Core i7-6700 Processor. Average completion times of over 10 minutes have not been recorded.

Figure \ref{fig:my-label} shows the time to compute the schema consequence for different schema sizes $\lvert S\rvert$ using \ensuremath{\texttt{critical}}\ and \ensuremath{\texttt{score}} . The parameters have been chosen to be small enough to accommodate for the high computational complexity of the \ensuremath{\texttt{critical}}\ approach.  This figure shows that \ensuremath{\texttt{score}}\ is orders of magnitude faster, especially on large schema sizes. The \ensuremath{\texttt{critical}}\ approach, instead, times out for schema sizes of over 33 triples.

Figure \ref{figure:scalabilityScore} shows the time to compute the schema consequence for different antecedent sizes $n_A$ and rule numbers $\lvert R\rvert$. The \ensuremath{\texttt{critical}}\ approach is not present in this figure, as it timed out in all the configurations. As this figure shows, the \ensuremath{\texttt{score}}\ approach can easily scale to a large set of rules.
Given the complexity of \texttt{SPARQL} query answering \cite{Perez2006Semantics}, we can also notice an exponential increase in computation time as more triples are added to the antecedent of a rule. In our experiment setup, the \ensuremath{\texttt{score}}\ approach scales to rules with antecedent sizes of up to 12 triples, before timing out.

\section{Related Work} \label{sec:relatedwork}

To the best of the authors' knowledge, our approach is the first to determine the applicability of inference rules to types of \texttt{RDF} triplestores specified by their schema, and to expand their schema with the potential consequences of such rules. Unlike related work on provenance paths for query inferences  \cite{GAD2017}, we do not explicitly model the dependencies between different rules. Instead, we compute their combined potential set of inferences by expanding the original schema on a rule-by-rule basis, through multiple iterations, following the basic principles of the chase algorithm.
We choose to follow a subset of the \texttt{RDF} data model, and not a simpler graph model such as \emph{generalised} \texttt{RDF} \cite{Hayes2014GeneralisedRDF}, to make our approach more applicable in practice, and compatible with existing tools. We pay particular attention to literals, as they are likely to feature prominently in IoT sensor measurements.

A possible application of our approach is to facilitate the sharing and reusing of inference rules. A relevant initiative in the IoT domain is Sensor-based Linked Open Rules (\texttt{S-LOR}) \cite{Gyrard2017linkedOpenRules}, which provides a comprehensive set of tools to deal with rules, including a rule discovery mechanism. By classifying rules according to sensor types, domain experts can discover and choose which inference rules are most relevant in a given scenario.
Our approach could automate parts of this process, by selecting rules applicable to the available data sources.
We refer to \cite{serrano2016review} for a comprehensive review of rule-based reasoning systems applicable to IoT.

Our approach to define a triplestore schema is related to a number of similar languages, and in particular to
Shape Expressions (ShEx) \cite{Prudhommeaux2014ShapeExpressions} and the Shapes Constraint Language (SHACL) \cite{2017SHACL}. The term \emph{shape}, in this case, refers to a particular constraint on the structure of an \texttt{RDF} graph.  ShEx and SHACL can be seen as very expressive schema languages, and computing schema consequences using such schemas would be impractical. In fact, each inference step would need to consider complex interpendencies between shapes and the values allowed in each triple element, and thus we would generate increasingly larger sets of contrainsts. The triplestore schema proposed in this paper is a simpler schema language and, if we disallow variables in the predicate position, can be modelled as a subset of both ShEx and SHACL.

\section{Conclusion}

As its main contribution, this paper presented two approaches to determine the applicability of a set rules with respect to a database schema (i.e.\ if the rule could ever match on any dataset modelled by the schema), by expanding such schema to model the potential facts that can be inferred using those rules, which we call \emph{schema consequence}.
This can be applied in IoT scenarios, where inference rules are used to aggregate sensor readings from diverse data sources in order to automate health and safety policies.  As the underlying data sources evolve, it is important to determine whether rules are still applicable, and what they can infer.

We focused on \texttt{RDF} triplestores, and on inference rules that can be modelled as \texttt{SPARQL} queries, such as \texttt{SPIN} and \texttt{SWRL}. To do so, we defined a notion of a \emph{triplestore schema} that constrains the type of triples allowed in a graph. This differs from the \texttt{RDF} \emph{schema} (\texttt{RDFS}) specification, which is designed to define vocabularies. While we provide an example on how to describe the schema of simple sensor networks, we think extending this approach to more expressive schema languages could be an interesting venue for future work.

The first of the two approaches that we presented is based on the existing notion of a critical instance; the second on query rewriting. We have theoretically demonstrated the functional equivalence of the approaches, as well as their soundness and completeness. Moreover, we have provided experimental evidence of the superior scalability of the second approach, which can be applied over large schemas and rulesets within seconds.

With this paper we intend to provide a first theoretical framework to reason about rule applicability. We hope that our approach will pave the way, on one hand, for efficient and meaningful policy reasoning in IoT systems, and on the other, for new and interesting rewriting-based schema reasoning approaches in the knowledge representation and databases research areas.

\bibliographystyle{splncs04}


\section{Introduction}
\label{intro}
The muscle-derived nicotinic acetylcholine receptor (nAChR) (PDB 2BG9) \citep{Unwin2005} is the most abundant neurotransmitter receptor at the neuromuscular junction (NMJ) in most vertebrates, including humans \citep{Albuquerque2009}. Within the postsynaptic membrane, nAChRs cluster in high densities (10,000 per $\mu$m$^{2}$) to properly activate the skeletal muscle \citep{Ramarao1998,Breckenridge1972}. As a major transmembrane protein, nAChR depends upon a highly specific lipid environment to maintain functionality. Lipids influence nAChR activity by affecting both function and organization. It is essential to understand how changes in lipid environment impact nAChR's structure and activity, given that lipids can change in response to aging and disease  \citep{Yadav2014}, and also vary across tissue and organism.

Over the past few decades, considerable progress has been made in uncovering lipid sensitivities associated with nAChR \citep{Criado1982}. In a majority of experiments, researchers have prioritized studying cholesterol over other membrane lipids. Early studies revealed that, when reconstituted into membrane mixtures, nAChR failed to conduct cations across the lipid bilayer unless cholesterol was present  \citep{Fong_Correlation_1986,Sunshine_Lipid_1992,Butler1993,Fong_Stabilization_1987,Corrie_Lipid_2002}.
More recently, researchers have examined the effects of membrane dynamics on nAChR organization \citep{Baenziger2015,Bruses2001,Marchand2002a,Oshikawa2003,Pato2008,Zhu2006a,Baenziger2017,Barrantes2007,Barrantes2000,Bermudez2010,Barrantes2010,Perillo2016,Wenz2005,Borroni2016,Unwin2017}. In-vitro studies \citep{Barrantes2007,Barrantes2010} indicate that nAChRs form larger aggregates upon cholesterol depletion. Experimental evidence suggest that cholesterol-rich lipid domains, known as lipid rafts, facilitate clustering of nAChRs \citep{Campagna2006,Marchand2002a,Pato2008}.
More specifically, after disrupting lipid raft formation, Zhu et al. observed a significant loss of nAChR clusters in-vitro \citep{Zhu2006a}. In the mature neuromuscular membrane, nAChRs are linked by the intracellular anchoring protein, rapsyn, which bridges receptors together at their bases  \citep{Zuber2013a}. According to fluorescent studies \citep{Marchand2002a}, lipid rafts mediate the association between rapsyn and neighboring nAChR molecules. In the mid-2000's, Willmann et al. and Stetzkowski-Marden et al. proposed that lipid rafts can stabilize receptor networks and may even provide a localized environment for nAChR \citep{Willmann2006,Stetzkowski-Marden2006}.

Domain formation occurs when a membrane is comprised of at least three lipid types: cholesterol, unsaturated fatty acids, and a molecule that interacts closely with cholesterol such as saturated fatty acids or sphingomyelin \citep{Feller_Acyl_2008,Yeagle2016115}.  The membrane separates into at least two domains, a liquid-ordered or ``raft'' phase containing cholesterol and saturated lipids/sphingomyelin, and a liquid-disordered phase containing unsaturated lipids. Polyunsaturated phospholipids make these domains more well-defined \citep{Levental_Polyunsaturated_2016}. Some results from atomistic simulations \citep{Sodt2014, Iyer2018} indicate that cholesterol may preferentially partition to the boundary between liquid-ordered phases and liquid-disordered phases composed of monounsaturated lipids.  In neuromuscular membranes, intrinsic domain formation is dependent upon several lipid species, including the widely influential omega-3 ($\omega$-3) fatty acids. One $\omega$-3 in particular, Docosahexaenoic acid (DHA), is prevalent in the native neuromuscular membrane and is strongly associated with flexible and well-defined domains  \citep{Turk2013,shaikh_dumaual_castillo_locascio_siddiqui_stillwell_wassall_2004}. Additionally, DHA is a major contributor to brain functioning, motor activity, and cardiac health; however, its specific effects on neuromuscular health are poorly understood \citep{12439486320170901,S000930840800032720080101,Georgieva2015}.

\begin{figure}[htp]
\centering
\includegraphics[width=.99\columnwidth,trim={0cm 2cm 0cm 0cm}]{nAChR2}
\caption[Nicotinic acetylcholine receptor (nAChR) structure]{{\bf Nicotinic acetylcholine receptor (nAChR) structure.} a) Atomistic and coarse-grained representations of Unwin et al's cryo-EM structure, PDB 2BG9. The nAChR is colored by subunit ($\alpha:green,\beta:purple,\delta:gray,\gamma:cyan$) and labeled by structure. The extracellular domain (ECD) is located above the bilayer and is critical for ligand-binding. The transmembrane domain (TMD) is positioned within the lipid bilayer, and the intracellular domain (ICD) is located in the cytoplasm. The coarse-grained model omits the ICD since it is poorly resolved and not necessary for this study. b) TMD from the extracellular perspective, with its four alpha helices labelled in a single $\alpha$ subunit. The outermost M4 helices closely interact with surrounding membrane lipids, while the M2 helices outline the ion pore; the M1 and M3 helices make up the body of the transmembrane domain \citep{Unwin2005}.}
  \label{fig:structure}
\end{figure}

Characterizing such complex lipid-protein interactions requires detailed information on each molecular structure. The nAChR is a member of a family of ion channels, known as pentameric ligand-gated ion channels (pLGICs), which has a number of recent structures \citep{Laverty2019,Laverty2017,Masiulis2019,Althoff2014,Hibbs2011,Morales2016,Baenziger2011,Corringer2012,S089662731630023X20160504,Prevost2012a,Sauguet2014}; however, only single structures have been determined, rather than dimers. pLGICs are composed of five subunits; each subunit is composed of an extracellular domain for ligand binding, a transmembrane domain (TMD), and an intracellular domain (ICD) (Figure 1A). The TMD is the embedded portion of pLGICs that interacts most often with surrounding membrane lipids; this region of pLGICs is composed of four alpha-helices, M1-M4, with the innermost, M2 helices forming the ion pore. The ICD, or cytoplasmic domain, is highly disordered, making it challenging to obtain information on its structure. Currently, there are two structures of nAChR pentamers resolved from x-ray crystallography (neuronal type) and cryo-electron microscopy (muscle-type) \citep{osti1330893, Unwin2005}.

Other closely related channels have been resolved with bound lipids. In 2017, the inhibitory neurotransmitter receptor, the gamma-aminobutyric acid (GABA$_A$) receptor, was crystallized with a cholesterol molecule protruding between its M1 and M3 helices \citep{Laverty2017}, similar to our previous prediction \citep{Henin2014}. Interestingly, x-ray structures of the glutamate-gated chloride channel (GluCl) revealed embedded POPC phospholipids in the same location \citep{Althoff2014}. Together, these data provide evidence for lipid-based modulation of pLGICs, which can potentially be applied to the gating of nAChR.

Structural information alone, however, is insufficient for answering questions about native cell membranes. For one, most structures are obtained under artificial conditions, using detergents or nanodisks with non-native lipid composition. Additionally, current structures do not represent oligomers, particularly in a liquid state. Fully atomistic molecular dynamics (MD) simulations \citep{brannigan, Cheng2009, Hnin_A_2014, Carswell_Role_2015} have complemented experimental investigations of pLGICs interacting with lipids, but they are limited in their ability to capture domain formation since atomistic lipids cannot diffuse over typical simulation microsecond time scales \citep{Ingolfsson2014,Bond2006,Parton2013,Goose2013,Scott2008}.  Furthermore, it is unfeasible to incorporate multiple nAChRs in simulations with atomistic resolution. Coarse-grained MD (CG-MD) simulations are run over longer length and time scales, making them suitable for exploring complex model membranes. CG-MD is widely applied in simulations of both lipid-protein binding and domain formation \citep{Bond2006,Scott2008,Parton2013,Goose2013}. Additionally, CG-MD can capture large-scale membrane phenomenon such as protein self-assembly and lipid-mediated oligomerization \citep{10.3389/fphys.2016.00494, BAADEN2013878}.

Recently, we \citep{Sharp2019} conducted CG-MD simulations of a single nAChR from the electric ray {\it Torpedo} in mixed membranes. Contrary to expectations, nAChR consistently preferred a local lipid environment rich in PUFAs rather than cholesterol, especially long-chained $\omega$-3s, such as DHA. While cholesterol occupied the transmembrane gaps of nAChR, PUFAs were even more likely to be embedded, regardless of their phospholipid headgroup.

The present study adopts a similar approach, with a particular focus on nAChR-associated clustering. Through molecular dynamics simulations, we investigate nAChR lipid preferences and clustering behavior in membranes with and without domains. For this study, we tested three major hypotheses: 1) Membrane organization affects nAChR boundary lipid specificity: when PUFA chains are prevented from forming PUFA rich domains, their prevalence among nAChR boundary lipids will be significantly reduced. 2) Domain formation will indirectly facilitate the clustering of nAChRs, by inducing partitioning preferences and restricting diffusion within the membrane. 3) Within a dimer, we will observe sequence preferences in facing subunits.
\section{Methods}
\subsection{System setup}

In order to isolate the role of lipid domain formation on nAChR-based lipid sorting and clustering, our simulations compared two different membranes with distinct phospholipid topology. Each membrane had 30\% cholesterol and 70\% phospholipids.  The phospholipids all had phosphatidylcholine (PC) headgroups, and half the total number of acyl chains were saturated (16:0) chains, while the other half were polyunsaturated (22:6) chains.  For homoacid membranes, all phospholipids had either two polyunsaturated acyl chains (didocosahexaenoylphosphatidylcholine or dDHA-PC) or two saturated chains (dipalmitoylphosphatidylcholine or DPPC), facilitating domain-formation. In the heteroacid membranes, all phospholipids were hybrid lipids with one polyunsaturated acyl chain and one saturated chain (1-palmitoyl- 2-docosahexaenoyl- phosphatidylcholine or PDPC), which are topologically unable to separate into PUFA-rich and PUFA-poor domains.

Lipids and proteins were modeled using the coarse-grained (CG) Martini 2.2 force field \citep{martini}.
Systems included between 1-4 nAChR molecules, derived from the Torpedo electric organ \citep{Unwin2005} (PDB 2BG9).  This structure is the only pLGIC structure obtained in a native membrane.

We converted protein structures into CG models using the Martini script "martinize.py", mapping four non-hydrogen atoms to one CG interaction. We constructed and assembled our protein-bilayer systems using  the Martini script, "insane.py", using a box sizes of 29x29x21 nm$^{3}$ for one-to-two proteins, 40x40x20 nm$^{3}$ for three proteins, and 44x44x22 nm$^{3}$ for four protein systems, respectively  \citep{Marrink2007}. Initially, the proteins were in a circle of about 13 nm. The receptors were evenly spaced, and their $\delta$ subunits were facing the same direction. Once simulations started, nAChRs shifted from their initial orientations. We ran 24 CG-MD simulations containing 1-4 nAChRs (3 replicas per system).

\subsection{Simulation details}
Simulations were run using the Martini 2.2 force field parameters and the Gromacs 5.1.2 simulation package, \citep{Marrink2007,Pronk2013} as in our previous work, \citep{Sharp2019}. Each simulation consisted of two steps: energy minimization and molecular dynamics. {For each system, we ran two consecutive energy minimizations for 10,000 steps.
Harmonic restraints between backbone atoms were imposed to preserve nAChR conformation.
More specifically, we applied an elastic force constant of 750 kJ/mol and set lower and upper bounds on the bond with using a bond length of 0.7 nm
 \citep{Marrink2007,Pronk2013}. The molecular dynamics simulations ran for 10-20 $\mu$s at a 0.025 ps time-step. Simulation temperature and pressure were kept constant at values of 323 K and a reference pressure of 1 bar. The isotropic pressure coupling compressibility constant was maintained at 3.0 $\times {10^{-5}}$ bar$^{-1}$.

\subsection{Analysis}

For a given nAChR, $n_{\mathrm{emb}}^{\alpha}$ is the total number of embedded lipids of lipid species $\alpha$, where embedded lipids satisfy the following criteria: the headgroup (in the case of cholesterol) or the terminal bead on the acyl chains (in the case of phospholipids) are within 10 {\AA}~of the M2 helices. Visual inspection indicated any lipids outside of this range were also outside of the TMD bundle.

Similarly, $n_{\mathrm{ann}}^{\alpha}$ is the total number of annular lipids of lipid species $\alpha$, where annular lipids satisfy the following criteria:  headgroup or the terminal bead on the acyl chains is between 10 {\AA}~and 35 {\AA}~of the M2 helices. This range of lipids corresponded to the third lipidation shell around nAChR.

(For phospholipids, each acyl chain was counted separately as a half-lipid, so it was possible for e.g. the sn-1 chain to be embedded and the sn-2 chain to be annular. )

Boundary lipid fractions for a given species $\alpha$ are defined as
\begin{eqnarray}
      f_{\mathrm{emb}}^{\alpha}&\equiv \frac{  n_{\mathrm{emb}}^{\alpha} }{n_{\mathrm{emb}} },\\
      f_{\mathrm{ann}}^{\alpha}&\equiv \frac{n_{\mathrm{ann}}^{\alpha}}{n_{\mathrm{ann}}}
    \label{eq:f}
  \end{eqnarray}
where $n_{\mathrm{emb}}$ and $n_{\mathrm{ann}}$ are the total number of embedded and annular lipids, respectively.

Two dimensional density distributions of boundary acyl chain species (B), $\tilde{\rho}_{B}(r_i,\theta_j)$ were calculated, as a function of radius, $r_i$, and angle $\theta_j$ projected onto the membrane.

  \begin{equation}

      \rho_{B}(r_i,\theta_j)= \frac{\left\langle n_{B}(r_i,\theta_j) \right\rangle}{r_i \Delta{r}\Delta{\theta}} \\

    \label{eq:R}
  \end{equation}

where $r_i = i \Delta r$ is the distance from the origin to the center of bin i, $\Delta r$ is the radial bin width, and $\Delta \theta$ is the angular bin width. $\left\langle n_{B}(r_i,\theta_j) \right\rangle$ is the time-averaged number of acyl chain beads of species B within bin i.

To quantify enrichment or depletion of a given acyl chain with respect to random distribution, the normalized density, $ \tilde{\rho}_{B}(r_i,\theta_j)$, was calculated.

  \begin{equation}

  \tilde{\rho}_{B}(r_i,\theta_j)=\frac{ \rho_{B}(r_i,\theta_j)}{x_{B}s_{B}~N_{L}/\langle L^{2}\rangle} \\

    \label{eq:Rt}
  \end{equation}

  where $s_B$ is the number of beads of lipid species B, $N_L$ is total number of lipids in a system, and $\langle L^{2}\rangle$ is the average projected box area, and $x_B$ is lipid B concentration. The expression does not take into consideration the protein footprint or undulations present within the system, and as such is an approximation.

The radial distribution function $g(r)$ of multiple nAChRs was calculated using the three dimensional distance $r$ between the centers of mass of the receptors for each of the 8000 frames, evenly distributed through the simulation.
It was derived from the distribution of pairwise separations, $P(r)$, by dividing by the expected separations in a random distribution : $g(r) = \frac{P(r)~dr}{r~dr}\times\int_{0}^{R} r dr$.
The height of the $g(r)$ peak corresponds to the amount of enrichment in probability for a given distance; for example, if g(10~nm) = 50, it is 50 times more probable that two monomers will be 10~nm apart than expected in a random distribution.

The total number of observed dimers $n_{d}$ is given by the number of pairs (across all analyzed proteins) where $r<100$ {\AA}. For each observed dimer, the closest subunit pair was determined, and the enrichment calculated as $F_{x,y}$:

\begin{equation}\label{eq:subunit}
F_{x,y} = \frac{{25}}{n_{d}}\left({n_{x,y}}+ n_{y,x}\right),
\end{equation}
where $n_{x,y}$ represents the number of dimers in which subunits $x$ and $y$ form the closest pair, and $1/25$ is the expectation for a random distribution

We visualized and imaged all simulation results using Visual Molecular Dynamics (VMD) \citep{HUMP96}.
\section{Results}
\begin{figure}[h]
\center
\includegraphics[width=\linewidth,trim={0cm 0cm 0cm 0cm}]{Fig1_cropped2.pdf}
\caption[Protein clustering in domain-forming (top row) and hybrid (bottom row) membranes.]{{\bf Protein clustering in domain-forming (top row) and hybrid (bottom row) membranes.} View of the membrane from the extracellular region at the final frame of 10 $\mu s$ simulations. nAChRs were colored by subunit ($\alpha:green,\beta:purple,\delta:gray,\gamma:cyan$), and lipids by acyl chain (DHA: white, saturated: blue, and cholesterol: red). Clustering of nAChRs on the borders of lipid rafts is visible in the domain-forming membranes.}
\label{fig:Figure4}
\end{figure}

\subsection{nAChR boundary lipid preferences}\label{sec:lipids}

\begin{table}[htbp]
  \centering
  \begin{tabular}{@{} c|cc|cc|c @{}}
     &\multicolumn{2}{c|}{mean embedded fraction}  &\multicolumn{2}{c|}{mean annular fraction} &bulk  \\
     & homoacid & heteroacid & homoacid & heteroacid  & \\
    \hline
    cholesterol & 0.29 (0.02)  & 0.42 (0.04) &0.23 (0.02)& 0.31 (0.01)& 0.30  \\

    PUFA  &0.61 (0.01)& 0.39 (0.04) & 0.52 (0.02) & 0.32 (0.01)&0.35\\

    saturated & 0.10 (0.05) & 0.19 (0.01) & 0.25 (0.02) & 0.37  (0.01)&0.35 \\

    \hline
  \end{tabular}
  \caption{Composition of nAChR Boundary Lipid chains in domain forming (homoacidic) and non-domain forming (heteroacidic) membranes.  Embedded and Annular chains are determined as described in Methods.  Averages are across systems, and across proteins in multi-protein systems.  Each protein is treated as a separate replica (n=30) and standard errors are shown in parentheses.  }
  \label{tab:bound}
\end{table}

In order to determine the significance of domain formation on enrichment of polyunsaturated acyl chains among boundary lipids, we calculated the fraction of embedded and annular boundary lipids in two types of membranes with equal amounts of cholesterol and saturated and polyunsaturated acyl chains.  In the domain-forming membranes, all phospholipids had either two polyunsaturated acyl chains or two saturated chains (homoacids), while in the non-domain forming membranes, all phospholipids had one polyunsaturated acyl chain and one saturated chain (heteroacids).

Distributions for $f_{emb}$ and $f_{ann}$ are shown in Figure \ref{fig:Fig2a} in heteroacids compositions without domains and in Figure \ref{fig:Fig2b} for homoacids compositions with domains, while mean values are given in Table \ref{tab:bound}.  For heteroacids, PUFA chains are slightly enriched among embedded chains and slightly depleted among annular chains, while the reverse trend is observed for saturated lipids.  Saturated chains are significantly depleted from the embedded lipids, but due to the topology of these lipids, a fully embedded polyunsaturated acyl chain must also contribute a saturated chain to the protein annulus. Distributions were remarkably consistent across simulations containing 1,2,3, or 4 proteins.

Figure \ref{fig:Fig2b} shows the results of the same analysis for homoacidic membranes.  Replacing heteroacids with homoacids substantially increases the fraction of polyunsaturated chains among both embedded and annular lipids, with the peak of the distribution shifting further to the right as more proteins are included (Figure \ref{fig:Fig2b}).

Together, these results are consistent with our previous \citep{Sharp2019} observation that polyunsaturated chains can displace cholesterol from embedded sites.  For heteroacids,  embedding a PUFA constrained the linked saturated chain to the nAChR annulus. In this case, cholesterol (which introduces no such constraint) was enriched among embedded lipids.  For homoacids, embedded PUFAs constrained a corresponding PUFA chain to the nAChR annulus, and here cholesterol was not enriched among embedded lipids.

The fraction of embedded chains that were polyunsaturated increased as more proteins were added to the homoacidic membrane. This could be consistent with each nAChR monomer in an oligomer blocking access to embedded lipid binding sites in the other monomers: flexible PUFA chains have multiple routes to access an embedded site, and rigid lipids have only one or two.

\begin{figure}[htp]
\centering
\includegraphics[width=\textwidth]{Boundary_a_ann.pdf}
\caption[nAChR boundary lipid preferences.]{{\bf Distribution of nAChR boundary cholesterol or acyl chain fractions in mixed membranes containing heteroacidic phospholipids.} Probability density distribution of fraction of embedded ($f_{emb}$) or annular lipids ($f_{ann}$) as defined in Eq. \ref{eq:f} are shown for 1 to 4 proteins.  Dashed lines represent expected boundary ratios for a randomly-mixed membrane, based on bulk lipid composition ($x_{CH} = 0.3$, $x_{SAT}=x_{PUFA}= 0.35$)}
  \label{fig:Fig2a}
\end{figure}

\begin{figure}[htp]
\centering
\includegraphics[width=\textwidth]{Boundary_b_ann.pdf}
\caption[nAChR boundary lipid preferences.]{{\bf Distribution of nAChR boundary cholesterol or acyl chain fractions in mixed membranes containing homoacidic phospholipids.} Probability density distribution of fraction of embedded ($f_{emb}$) or annular lipids ($f_{ann}$) as defined in Eq. \ref{eq:f} are shown for 1 to 4 proteins.  Dashed lines represent expected boundary ratios for a randomly-mixed membrane, based on bulk lipid composition ($x_{CH} = 0.3$, $x_{SAT}=x_{PUFA}= 0.35$).} \label{fig:Fig2b}
\end{figure}

\begin{figure}[htp]
\centering
\includegraphics[width=\textwidth]{Polar_plot.pdf}
\caption[Polar plot.]{{\bf Enrichment or depletion of lipid species surrounding a single nAChR.} Heatmaps represent normalized lipid densities ($\tilde{\rho}_{B}(r_i,\theta_j)$) as defined in Eq. \ref{eq:Rt}, on a scale of -3 (blue, indicating depletion) to 1 (red, indicating enrichment). Densities were averaged over the last 8 $\mu s$ of the 10 $\mu s$ simulations and across three replicas.} \label{fig:Fig2c}
\end{figure}

The two dimensional density distribution of each lipid species from the center of a single nAChR is shown in Figure \ref{fig:Fig2c}. For all three lipids, there was a five-fold symmetry of densities surrounding nAChR, with lipid preferences determined by transmembrane helix, rather than by subunit. In heteroacid mixtures, cholesterol is enriched near the subunit interfaces, and even buried more deeply within the $\gamma/\alpha$ interface, with saturated acyl chains packed just outside cholesterol.  The outermost M4 helices, however, were packed with PUFAs; PUFA chains also diffused throughout the TMD bundle.  In homoacid membranes, saturated lipids and cholesterol were depleted among embedded lipids, but they maintained the same helix association observed in heteroacid membranes. PUFAs were especially enriched in domain-forming membranes, with high densities around all four transmembrane helices.

\subsection{nAChR clustering in the presence and absence of lipid domains} \label{sec:clustering}
Here, we repeatedly observed spontaneous formation of receptor dimers in the simulations containing multiple proteins (Figure \ref{fig:Figure4}).
To investigate the role of lipid domain formation in driving nAChR oligomerization, we calculated the radial distribution function for the pairwise distances between centers of mass, as shown in Figure \ref{fig:Figure6}. For both heteroacid and homoacid membranes, we observed a peak at $r\sim7.5$~nm, corresponding to dimerization.    The differences in profiles between domain and non-domain forming membranes were primarily quantitative: most significantly, the peak corresponding to dimers was substantially higher for 2 and 3 proteins in homoacids than in heteroacids. For two proteins, the distribution for homoacids was shifted to the left, relative to the distribution for heteroacids, indicating that the peak for dimerization was at a shorter distance in domain-forming membranes. This difference is consistent with results indicating a role for domain formation in aggregation and clustering of nAChRs \citep{Barrantes2007,Oshikawa2003,Pato2008}. The simplest explanation of this difference is the higher effective concentration of proteins when domains are present: all nAChRs are corralled in a single liquid-disordered domain, with approximately half the area of the overall membrane.

\begin{figure}[htp]
\centering
\includegraphics[width=1.0\linewidth]{Fig3_gr_new_bins}
\caption[Pairwise distance distribution across multi-protein systems.] {{\bf Radial distribution function g(r)  of pairwise protein center-of-mass distances, across multi-protein systems.}  Data collection began at 2 $\mu s$ into each trajectory and ended at 10 $\mu s$, respectively. The peak between 7-10 nm corresponds to dimer formation.  }
\label{fig:Figure6}
\end{figure}

\subsection{Closest subunits across dimerizing proteins}
In order to determine whether nAChR dimers were more likely to form with specific subunit interactions, we determined the closest interacting subunit pair for nAChR dimers within each frame. To ensure that dimers, rather than larger oligomers, were being analyzed, we excluded trimers and tetramers from the analysis and only considered two protein systems, with increased sampling. Figure ~\ref{fig:Figure7} shows the amount of enrichment for each possible subunit pairing, relative to an expected random distribution. Results were very sensitive to the use of domain-forming compositions. In heteroacidic membranes, the $\alpha_{\delta}$ subunit formed a monomer-monomer interaction with the $\beta$ subunit, while the $\alpha_{\gamma}$ most favorably interacted with the $\delta$ subunit. In domain-forming membranes,  the $\delta-\alpha_{\delta}-\gamma$ interface was two to five times as likely to pair with a $\alpha_{\delta}$ subunit compared to a random distribution, but substantially more dispersion was detected than in non-domain forming membranes. This difference suggests that nAChRs have preferred dimerization orientations, but membrane organization can reorient receptors. Although $\delta$ subunits are linked by a disulfide bond at the NMJ \citep{Chang1977}, our simulations suggest that without this link, $\delta$ subunits do not form the closest pair among dimerizing proteins.

\begin{figure}[htp]
	\center
	\includegraphics[height=3in]{Fig4_new}
\caption[Subunit pairs among dimerizing proteins]{{\bf Most probable subunit pairs among dimerizing proteins.} Heat-map showing the likelihood that two subunits were closest together among dimerizing proteins (100 {\AA} threshold). Distinct subunits are each expected to form the closest pair 8 \% of the time, while identical subunits are each expected to dimerize only 4 \% of the time. Each color on the heatmap represents multiplicative depletion or enrichment relative to the pair's expectation value, ranging from complete depletion (red) to no enrichment (white) to five-fold enrichment (blue).  }
\label{fig:Figure7}
\end{figure}

\section{Discussion}

The nAChR is one of the most well-studied, fundamental pLGICs for understanding human cognition, memory, and muscle contraction \citep{Gotti1997}. As an integral membrane protein, the function and organization of nAChR is strongly dictated by its surrounding lipid environment. DHA is an $\omega$-3 polyunsaturated fatty acid abundant in synaptic membranes and the Torpedo electric organs, which are both native nAChR membranes. With its six double bonds, DHA is considered highly disordered and can induce domain formation in membranes \citep{S000930840800032720080101}. We previously observed \citep{Sharp2019} favorable interactions between a single nAChR and DHA-rich, cholesterol-poor domains using coarse-grained simulations. Here, we have extended this approach to investigate the role of lipid topology and domain-formation on boundary lipids of individual nAChRs. We have also conducted the first simulations containing multiple nAChRs, which has allowed us to observe spontaneous dimerization.

While DHA is implicated in human health and disease, \citep{12439486320170901} experimental studies considering its interactions with nAChR have exclusively considered its free-fatty acid form (FFA), \citep{Antollini2016}
rather than as an acyl chain component of a phospholipid. Application of $\omega-3$ FFAs causes a significant reduction in open times observed through single-channel recordings  \citep{Bouzat1993}. Here, we observe a substantial effect of lipid topology on both embedded and annular lipids: DHA chains in homoacidic phospholipids are far more likely to be found as either annular or embedded boundary lipids.  We previously observed only quantitative effects of swapping PE with PC on partitioning, but the headgroup does serve to anchor the lipid at the membrane/protein interface  \citep{Sharp2019}.  DHA in its FFA form (without a headgroup) may diffuse into an open nAChR pore, blocking the channel.

We find that, consistent with coarse-grained MD simulations using one nAChR \citep{Sharp2019}, multiple nAChRs continue to prefer the liquid-disordered phase containing long-chain $\omega$-3 fatty acids. While the number of nAChRs in the system did not affect the partitioning profile in these simulations, it did affect the composition of embedded lipids. Our results are consistent with each nAChR monomer in an oligomer blocking access to embedded lipid binding sites in the other monomers, suggesting an intriguing coupling between specific binding and membrane organization.

Interestingly, upon removing membrane organization, embedded lipids cluster around specific transmembrane helices in a five-fold symmetry around nAChR. This finding suggests that intrinsic lipid preferences are primarily helix dependent, rather than subunit dependent. In homoacid membranes, a lipid preference for DHA was observed across all transmembrane helices, with shells of PUFAs found even at the border of the liquid-ordered phase. Although saturated acyl chains and cholesterol were generally depleted around nAChR, the highest densities for both lipids were found around the M1 and M3 helices, as seen in heteroacid membranes.

In native membranes, nAChR dimers can be stabilized by a disulfide bond between $\delta$ subunits \citep{Chang1977}.  An early controversy \citep{Anholt1980,Ruechel1981,Zingsheim1982,Schindler1984} concerned whether the disulfide bond was necessary for dimer formation. There is no mechanism for covalent bonds between monomers in these coarse-grained simulations. All dimers were stabilized by non-covalent interactions, consistent with the results of \citep{Ruechel1981,Schindler1984}. We observed far more stable dimers in the homoacid mixtures than the heteroacid mixtures, which would be consistent with high sensitivity to experimental conditions. We did not observe a consistent $\delta-\delta$ preference for interfacing subunits in either homoacid or heteroacid mixtures. Schindler {\it et al} (1984) observed an apparent gain-of-function for single-channels within dimers relative to monomers, regardless of disulfide linking \citep{Schindler1984}, suggesting that lipid modulation of oligomerization could also provide a pathway for modulating single channel function.

\begin{acknowledgements}
GB was supported by research grants NSF MCB1330728 and NIH P01GM55876-14A1. GB and LM were also supported through a grant from the Research Corporation for Scientific Advancement. This project was supported with computational resources from the National Science Foundation XSEDE program through allocation NSF-MCB110149, a local cluster funded by NSF-DBI1126052, the Rutgers University Office of Advanced Research Computing (OARC) and the Rutgers Discovery Informatics Institute (RDI2), which is supported by Rutgers and the State of New Jersey. We are grateful to Dr. J\'{e}r\^{o}me H\'{e}nin for his helpful suggestions throughout this study.
\end{acknowledgements}

\section{Compliance with Ethical Standards}
{\bf Conflict of interest:} The authors declare that they have no conflict of interest. This research was supported in part by the National Science Foundation, the National Institutes of Health, and the Research Corporation for Scientific Advancement.

\bibliographystyle{spbasic}