11institutetext: Lun Ai 22institutetext: Department of Computing, Imperial
College London, London, UK
22email<EMAIL_ADDRESS>33institutetext: Stephen H. Muggleton
44institutetext: Department of Computing, Imperial College London, London, UK
44email<EMAIL_ADDRESS>55institutetext: Shi-Shun Liang
66institutetext: Department of Life Science, Imperial College London, London,
UK
66email<EMAIL_ADDRESS>77institutetext: Geoff S.
Baldwin88institutetext: Department of Life Science, Imperial College London,
London, UK
88email<EMAIL_ADDRESS>
# Boolean matrix logic programming for active learning of gene functions in
genome-scale metabolic network models

Lun Ai  Stephen H. Muggleton  Shi-Shun Liang  Geoff S. Baldwin

###### Abstract

Techniques to autonomously drive research have been prominent in Computational
Scientific Discovery, while Synthetic Biology is a field of science that
focuses on designing and constructing new biological systems for useful
purposes. Here we seek to apply logic-based machine learning techniques to
facilitate cellular engineering and drive biological discovery. Comprehensive
databases of metabolic processes called genome-scale metabolic network models
(GEMs) are often used to evaluate cellular engineering strategies to optimise
target compound production. However, predicted host behaviours are not always
correctly described by GEMs, often due to errors in the models. The task of
learning the intricate genetic interactions within GEMs presents computational
and empirical challenges. To address these, we describe a novel approach
called Boolean Matrix Logic Programming (BMLP) by leveraging boolean matrices
to evaluate large logic programs. We introduce a new system, $BMLP_{active}$,
which efficiently explores the genomic hypothesis space by guiding informative
experimentation through active learning. In contrast to sub-symbolic methods,
$BMLP_{active}$ encodes a state-of-the-art GEM of a widely accepted bacterial
host in an interpretable and logical representation using datalog logic
programs. Notably, $BMLP_{active}$ can successfully learn the interaction
between a gene pair with fewer training examples than random experimentation,
overcoming the increase in experimental design space. $BMLP_{active}$ enables
rapid optimisation of metabolic models to reliably engineer biological systems
for producing useful compounds. It offers a realistic approach to creating a
self-driving lab for microbial engineering.

###### Keywords:

Computational Scientific Discovery; Synthetic Biology; Active Learning;
Inductive Logic Programming; Matrix.

††journal: Machine learning

##  1 Introduction

A priority in scientific research is the development of autonomous systems to
enhance the scientific discovery process. Currently, most efforts in this
domain focus on either laboratory automation or machine learning applied to
scientific datasets. Sub-symbolic approaches usually demand a large number of
annotated data and extensive trial-and-error. It is often difficult to modify
the structure of these machine learning systems or extract human-interpretable
information. The application of these systems in an experimental workflow is
limited by data efficiency, flexibility and interpretability. There is a
resurgence of interest in incorporating AI in experimental automation and
scientific understanding simultaneously. In this work, we investigate an
integration of techniques for planning experiments, optimising resource
utilisation and interpretable machine learning to support biological
discovery.

In the pursuit of understanding natural phenomena, scientists often formulate
plausible theories and design experiments to test competing hypotheses. This
process typically involves a careful selection of experimental conditions and
the logical inference mechanism of abduction. Experimentation can be optimised
through active learning, where computational methods are employed to
strategically select informative training examples to combat resource
limitations. On the other hand, Inductive Logic Programming (ILP) ILP1991
offers an automated approach to abduction, representing observations,
hypotheses and background knowledge through interpretable logic programs. In
ILP, hypotheses are learned to explain observational data with respect to the
background knowledge. “Askable” hypotheses, known as abducibles, that might
explain contradictions between observations and known background knowledge are
first assumed to hold ase_progol . The resulting extended knowledge, enriched
with abduced hypotheses, is subsequently verified against observations TCIE .

Our study explores the applicability of abductive reasoning and active
learning to comprehensive biological systems. Given that biological
relationships are commonly described logically, ILP is particularly adept at
operating on biological knowledge bases. The integration of abductive
reasoning and active learning via ILP was successfully demonstrated in the
context of biological discovery by the prominent Robot Scientist
King04:RobotScientist . The Robot Scientist performed active learning by
strategically selecting key experiments to achieve more data- and cost-
effective gene function learning than random experimentation. However, this
demonstration was limited to only 17 genes in the aromatic amino acid pathway
of yeast. In the present study, we advance the scope of logic-based gene
function learning by looking at Escherichia coli (E. coli). The genome-scale
metabolic model (GEM), iML1515 iML1515 , encompasses 1515 genes and 2719
metabolic reactions of the E. coli strain K-12 MG1655, a versatile host
organism for metabolic engineering to produce specific compounds. Notably, the
iML1515 model represents a 90-fold increase in complexity compared to the
aromatic amino acid pathway investigated by the Robot Scientist.

Despite the extensive mapping of the E. coli genome to metabolic functions in
iML1515, it was shown to yield inaccurate phenotype predictions
bernstein_evaluating_2023 . A classical method to address errors in metabolic
pathways involves auxotrophic growth experiments, where specific genes are
deleted to render the organism incapable of synthesizing essential compounds
beadle_genetic_1941 . For accurate simulation of auxotrophic mutant
experiments, the computational representation of GEMs must faithfully depict
the underlying network structure and effectively facilitate the exploration of
metabolic pathways.

Considering these requirements, Petri nets emerge as suitable modelling tools
for GEMs. Petri nets are a class of directed bipartite graphs with diverse
applications in Computer Science and Biology reisig_understanding_2013 ;
sahu_advances_2021 . Petri nets can represent a structured knowledge base
while capturing dynamic behaviours in a system. In contrast, logic programs
allow reasoning over relationships between entities but cannot directly
reflect the changes in concepts and entity states. To integrate logic
reasoning with Petri nets, we explore a more holistic approach to learning
metabolic networks that combines recursive logic programs with pathway
reachability analysis.

Figure 1: The hypotheses and instances are represented as logic programs to
describe gene functions (gene-enzyme associations) and auxotrophic mutant
experiment conditions respectively. $BMLP_{active}$ encodes the GEM iML1515
and simulates the auxotrophic mutant experiment outcomes. It actively consults
a data source to request ground truth labels that minimise the expected
experimental cost based on a user-provided cost function when pruning
hypotheses. The underlying BMLP iteratively refutes gene function hypotheses
inconsistent with labelled training examples. While $BMLP_{active}$ to date
only learns from synthetic data, a laboratory robot can be integrated to
perform selected experiments.

Learning complex biological systems like GEMs presents significant challenges
in terms of computational complexity and empirical feasibility. To address
these challenges, we present a new system $BMLP_{active}$ (Figure 1) based on
a novel approach called Boolean Matrix Logic Programming (BMLP) and an
implementation in SWI-Prolog. We use logic programs to encode the biochemical,
genetic, and genomic relationships in iML1515. While simulating auxotrophic
mutant experiments, logical inferences to predict phenotypic effects are
accelerated by boolean matrices via BMLP. This speeds up the runtime for
predicting phenotypic effects by a factor of 170 compared to SWI-Prolog
without BMLP on a single computation thread. We implement active learning to
select auxotrophic mutant experiments that minimise the expected experimental
cost via a user-defined cost function. Empirically, we delete gene functions
in iML1515 and use $BMLP_{active}$ to learn hypotheses to recover the deleted
gene functions. $BMLP_{active}$ successfully restores these removed gene
functions through active learning while reducing 90% the optional nutrient
substance cost required by randomly selected experiments. It additionally
lowers the number of experimental data needed to recover these gene functions
compared to random experiment selection, reducing the sample complexity. When
operating within a finite budget, our approach can deliver optimal
experimental outcomes, whereas studies employing random experiment selection
might not reach completion.

Furthermore, we apply $BMLP_{active}$ to learning digenic functions, which
holds significant implications for drug development and therapeutic
interventions costanzo_global_2016 . These digenic interactions are related to
isoenzymes - two enzymes responsible for the same reaction. The intricate
relationship between complex genotypes and phenotypes remains largely
unexplored in most organisms costanzo_global_2019 . Additionally, some digenic
interactions are dynamic and can vary depending on the growth conditions
bernstein_evaluating_2023 . A comprehensive understanding of digenic
interactions would expedite strain engineering efforts. We show that
$BMLP_{active}$ can converge to the correct gene-isoenzyme mapping with as
little as 20 training examples. This shows a significant promise for
$BMLP_{active}$ to address complex genetic interactions on the whole genome
level.

Our contributions are:

  1. 1.

We describe a novel problem of Boolean Matrix Logic Programming (BMLP).

  2. 2.

We implement a BMLP algorithm, BMLP-IE, in SWI-Prolog to predict phenotypic
effects by computing reachability queries in GEMs based on a subclass of Petri
net, which empirically improves simulation runtime by 170-fold compared to
SWI-Prolog without BMLP-IE.

  3. 3.

We built a BMLP system, $BMLP_{active}$, which uses BMLP-IE to actively select
experiment instances for learning gene functions in a GEM.

  4. 4.

We present a theoretical sample complexity ratio between active learning and
random experiment sampling and empirically highlight the effect of actively
selecting examples in reducing the total training examples needed.

  5. 5.

We empirically show that $BMLP_{active}$ reduces 90% optional nutrient
substance cost and the number of experiments needed for learning gene
functions compared to randomly sampled experiments.

  6. 6.

Our empirical results demonstrate that $BMLP_{active}$ can be devised to learn
gene-isoenzyme associations in GEM with a faster convergence to the correct
gene-isoenzyme mapping than randomly sampled experiments.

The rest of this paper is organised as follows. In Section 2, we summarise
work on machine learning of genome-scale metabolic networks, matrix algebra,
Petri nets and active learning. In Section 3, we introduce key concepts in
matrix algebra and an elementary class of Petri nets that we use to model a
GEM, called one-bounded elementary nets (OENs). In Section 4, we describe a
transformation from an OEN to a subclass of recursive logic programs, linear
and immediately recursive datalog. In Section 5, we introduce our BMLP
framework to evaluate this class of logic programs. In Section 6, we describe
our active learning approach based on prior work ase_progol ; progol . In
Section 7, we explain the $BMLP_{active}$ implementation for encoding iML1515,
GEM simulation and active learning gene functions. Section 8 includes our
empirical results on simulation runtime, phenotype prediction, active learning
gene functions and gene-isoenzyme mapping. In Section 9, we explain the
current limitations of $BMLP_{active}$ and discuss future work.

##  2 Related work

###  2.1 Active learning in computational scientific discovery

Active learning techniques provide machine learning algorithms with more
control over the training examples they use cohn_improving_1994 . A typical
approach to active learning is via membership queries where the machine
learner queries a data source to label instances according to some strategy.
An active learner can be implemented to explore the version space which
contains all hypotheses that are expressible by a hypothesis language and are
consistent with a set of training examples mitchell_generalization_1982 .
Instances are chosen according to a measure of informativeness such as entropy
shannon_mathematical_1963 , diameter tosh_diameter_based_2017 , the size
mitchell_generalization_1982 ; dasgupta_coarse_2005 and the shrinkage
hocquette_how_2018 of the version space for eliminating competing hypotheses.
For a binary classification task, when learners are allowed to look at $n$
unlabelled instances to compare, an active learner has $n$ times larger
probability of selecting an instance with maximal entropy than a passive
learner hocquette_how_2018 . The general principle for a binary active learner
is to label instances closest to being consistent with half of the hypotheses
in the version space mitchell_generalization_1982 . These instances reject up
to half of the hypothesis space regardless of their classifications. In the
optimal case when all instances are maximally discriminative, the target
hypothesis can be found via a binary search with depth logarithmic of the
hypothesis space size angluin_queries_2001 . $BMLP_{active}$ makes membership
queries but our objective is to not only reduce the number of labels but also
minimise the expected cost of experimentation.

Computational systems such as the BACON systems langley_scientific_1987 , the
LAGRANGE systems todorovski_integrating_2006 , recent equation discovery
systems brunton_discovering_2016 ; guimera_bayesian_2020 ; petersen_deep_2020
are primarily designed to formulate symbolic hypotheses from experimental
results. However, one cannot directly employ these systems for experimental
planning. The automation of experimental design was investigated by the
pioneering Robot Scientist King04:RobotScientist as an endeavour to bring
together logical reasoning, active learning and laboratory automation. The
Robot Scientist automatically proposes explanatory hypotheses, actively
devises validation experiments, performs experiments using laboratory
robotics, and interprets empirical observations. Abduction and active learning
are supported by ASE-Progol ase_progol based on the ILP system Progol progol .
Experiments were actively selected to minimise the expected cost of
experimentation for learning gene functions in the aromatic amino acid
synthesis pathway of yeast. It was combined with a laboratory robot to
automatically reconstruct the yeast metabolic pathway model and the active
learning strategy significantly outperformed random experiment selection
King04:RobotScientist . Later, a high-throughput laboratory automation
platform, Adam king_automation_2009 extended this approach by constructing a
logical model of yeast consisting of 1200 genes. In comparison, we apply
$BMLP_{active}$ to a genome-scale metabolic network, which is significantly
larger than the models in King04:RobotScientist ; king_automation_2009 .

Recently, large language models (LLMs) have been explored to automate
experimental design via user-provided prompt texts boiko_autonomous_2023 ;
bran_augmenting_2023 . LLMs are generative models of natural languages based
on neural networks that are trained from textual data. In
boiko_autonomous_2023 ; bran_augmenting_2023 , LLMs served as experiment
planners and coding assistants. These platforms use LLMs for planning
experimental steps using publicly available knowledge bases, writing and
interpretation programs to control hardware in a laboratory environment and
analysing collected experimental data. Their advantage is accessing rich
information online and interacting with users while designing experiments.
However, the consistency and faithfulness of their outputs cannot be
guaranteed. It has been stressed that autonomous experimentation by AI should
be employed with caution without the means to identify the source of error and
explain observations ren_autonomous_2023 . While $BMLP_{active}$ cannot yet
interact with biologists in natural language, interpreting and verifying its
output are straightforward due to its logical representation and computation.

###  2.2 Learning genome-scale metabolic network models

GEMs contain known metabolic reactions in organisms, serving as functional
knowledge bases for their metabolic processes. There has been an increasing
interest in integrating machine learning methods with GEMs to improve the
predictions of gene-phenotype correlations within biological systems
angione_human_2019 ; sen_deep_2021 . This often requires constraint-based
modelling techniques for simulating GEMs rana_recent_2020 . Flux balance
analysis (FBA) palsson_systems_2015 is the most common mathematical approach
to solve flux distributions in genome-scale models under the steady-state
condition. FBA uses linear programming to model fluxes in a chemical reaction
network represented by a matrix that contains the stoichiometric coefficients
of metabolites in reactions. Ensemble learning frameworks such as
wu_rapid_2016 ; oyetunde_machine_2019 ; heckmann_machine_2018 obtain support
vector machines, decision trees and artificial neural networks to represent
hidden constraints between genetic factors and metabolic fluxes. These systems
cannot identify inconsistencies between simulation results and experimental
data to self-initialise constraints to aid learning rana_recent_2020 . On the
other hand, mechanistic information that reflects whole-cell dynamics has been
incorporated to address these concerns, for instance in artificial neural
networks to tune their parameters yang_white_box_2019 . The hybrid approach in
faure_neural-mechanistic_2023 embeds FBA within artificial neural networks
based on custom loss functions surrogating the FBA constraints. Recent
research also explores autoencoders for learning the underlying relationship
between gene expression data and metabolic fluxes via fluxes estimated by FBA
sahu_advances_2021 . Currently, machine learning of genome-scale metabolic
networks is significantly limited by the availability of training data
sen_deep_2021 . In contrast, $BMLP_{active}$ does not use FBA so it is not
dependent on the metabolic flux constraints. $BMLP_{active}$ also performs
active learning to reduce the cost and the number of experiments. This is a
notable advantage considering finite experiment resources and insufficient
experimental data.

###  2.3 Matrix algebraic approaches to logic programming

In our BMLP framework, we use boolean matrices to evaluate fixpoints of large
recursive datalog programs. Obtaining the least Herbrand model of recursive
datalog programs can be reduced to computing the transitive closure of boolean
matrices peirce_collected_1932 ; copilowish_matrix_1948 . Fischer and Meyer
fischer_boolean_1971 study a divide-and-conquer boolean matrix computation
technique by viewing relational databases as graphs. A similar approach is
explored by Ioannidis ioannidis_computation_1986 for computing the fixpoint of
recursive Horn clauses. An ILP system called DeepLog
muggleton_hypothesizing_2023 employs boolean matrices for learning recursive
datalog programs. It repeatedly computes square boolean matrices to quickly
choose optimal background knowledge for deriving target hypotheses. The
mapping from biological knowledge bases to propositional logic programs is
explored in sato_boolean_2021 ; sato_differentiable_2023 where logic programs
are represented as boolean matrices by considering gene regulatory networks as
boolean networks. In contrast, we encode a GEM as a Petri net and a recursive
datalog program for computing reachability queries.

Recent work primarily studies the mapping of logic programs to linear
equations over tensor spaces. SAT problems are investigated in
lin_satisability_2013 under a linear algebra framework. Truth values of domain
entities, logical relations and operators for first-order logic can be
evaluated using tensors-based calculus grefenstette_towards_2013 . Computing
algebraic solutions can approximate recursive relations in first-order Datalog
sato_linear_2017 . Based on this approach, abduction can be performed by
encoding a subset of recursive datalog programs in tensor space
sato_abducing_2018 . Tensor representations allow differentiable
infrastructures such as neural networks to perform probabilistic logical
inferences cohen_tensorlog_2020 . This can be done using dynamic programming
to propagate beliefs for stochastic first-order logic programs. The
probabilistic parameters can be optimised using deep learning frameworks
yang_differentiable_2017 . Limitations have been identified for solving
certain recursive programs due to the difficulty of computing some arithmetic
solutions whereas these can be solved by iterative bottom-up evaluations
sato_linear_2017 . In addition, algebraic approximations are usually
calculated in linear algebraic methods where the correctness of the solutions
cannot be guaranteed. However, our BMLP approach does not have these issues
since it uses boolean matrices for iterative bottom-up computation.

###  2.4 Petri nets in logic programming

Little attention has been paid to representing and simulating Petri nets with
logic programs. With logic programming, techniques to execute, analyse and
learn Petri nets can be implemented in a single programming language
domenici_petri_1990 . The early work domenici_petri_1990 describes places and
transitions of a Petri net in terms of n-ary ground facts. Srinivasan et al.
srinivasan_knowledge_2012 attempt to learn n-ary ground facts that represent
Petri net transitions in the ILP ILP1991 system Aleph aleph . Few studies have
looked at encoding Petri nets by Answer Set Programming (ASP)
gelfond_stable_1988 . Behrens and Dix behrens_model_2007 use Petri nets for
analysing interactions in multi-agent systems with ASP programs. Anwar et al.
anwar_encoding_2013 provide an ASP encoding of Petri nets with coloured tokens
to model the behaviours of biological systems. This encoding is adapted by
Dimopoulos et al. dimopoulos_encoding_2020 for simulating reversible Petri
nets that allow transitions to fire in both directions. These systems all use
high-level symbolic representations of Petri nets as deductive knowledge
bases. Logic programming approaches rely on navigating background knowledge to
determine the consistency of hypotheses. However, non-high-throughput
deductive inferences limit these systems’ applicability for dealing with
increasingly complex knowledge bases such as GEMs. $BMLP_{active}$ addresses
this limitation by utilising a boolean matrix encoding of GEMs and boolean
matrix operations to accelerate inferences.

##  3 Preliminaries

###  3.1 Matrix algebra for datalog

We introduce datalog programs and concepts relating to our BMLP framework
based on LP_book ; ILP_foundation . A variable is a character string beginning
with an uppercase letter. A function, predicate or constant symbol is a
character string starting with a lowercase letter. A term is either a
variable, a constant symbol, or a function symbol attached to a tuple of
terms. A ground term contains no variables. An atom is a predicate symbol
applied to a tuple of terms e.g. $r(a,Y)$ where $r$ is a predicate symbol, $a$
is a constant symbol and $Y$ is a variable. An atom is ground when it contains
only ground terms. A Horn clause is a disjunctive clause with at most one
positive literal. A definite clause is a Horn clause with exactly one positive
literal and a definite program is a set of definite clauses. A query is a
user-posted Horn clause without any positive literal. A general logic program
engine such as a Prolog interpreter can evaluate a query by carrying out a
refutation proof of its negation via SLD resolutions (Selective Linear
Definite clause resolution). A unit ground clause (fact) is a definite clause
with no body literal and no variables. Datalog programs are definite programs
without function symbols. A first-order definite clause has the form,
$A_{0}\leftarrow A_{1},...,A_{n}$ where $A_{i}=r_{i}(x_{1},...,x_{k})$,
$r_{i}$ is a predicate symbol, $x_{j}$ is a term. The arity of a function or
predicate symbol is its number of arguments, e.g. $r_{i}/k$ indicates $r_{i}$
has arity $k$. The literal $A_{0}$ is called the head of the clause and the
literals $A_{i}$ are the body.

The Herbrand base $\mathcal{B_{P}}$ is the set of all unit ground clauses
provable from the predicate and constant symbols in a first-order definite
program $\mathcal{P}$. The Immediate Consequence Operator $\mathcal{T_{P}}$
van_emden_semantics_1976 is a mapping over subsets of $\mathcal{B_{P}}$
defined as:

$\displaystyle\forall
I\subseteq\mathcal{B_{P}},\mathcal{T_{P}}(I)=\\{\alpha\in\mathcal{B_{P}}\,|\,\alpha\leftarrow
B_{1},...,B_{m}\,(m\geq 0)$ $\displaystyle\text{ is a ground instance of a
clause in }\mathcal{P}\text{ and }\\{B_{1},...,B_{m}\\}\subseteq I\\}.$

$\mathcal{T_{P}}$’s least fixpoint is $\mathcal{T_{P}}{\uparrow^{\omega}}$ =
$\cup_{n\geq 0}\mathcal{T_{P}}{\uparrow^{n}}$ where
$\mathcal{T_{P}}{\uparrow^{0}}=\emptyset$ and
$\mathcal{T_{P}}{\uparrow^{n}}=\mathcal{T_{P}}(\mathcal{T_{P}}{\uparrow^{n-1}})$.
A clause entailed by $\mathcal{P}$ is written as $\mathcal{P}\models\alpha$.
The minimal Herbrand model $M(\mathcal{P})$ contains all unit ground clauses
that are entailed by $\mathcal{P}$ and $M(\mathcal{P})$ =
$\mathcal{T_{P}}{\uparrow^{\omega}}$.

A recursive program has a body literal that appears in the head of a clause.
Specifically, we focus on linear and immediately recursive datalog as simple
recursive datalog programs. A linear and immediately recursive datalog program
has a single recursive clause with the recursive predicate appearing in the
body only once ioannidis_computation_1986 . This subset of recursive datalog
programs can represent the relational structure of reaction pathways. This is
further elaborated in Section 4. An example of a linear and immediately
recursive datalog program involving one non-recursive predicate is:

$\displaystyle r_{2}(X,Y)\leftarrow r_{1}(X,Y).$ $\displaystyle
r_{2}(X,Y)\leftarrow r_{1}(X,Z),r_{2}(Z,Y).$

A recursive datalog program $\mathcal{P}$ can be written as matrices for
bottom-up evaluation ceri_datalog_1989 . When $\mathcal{P}$ is linear and
immediately recursive, it can be described by linear equations
sato_linear_2017 :

$\displaystyle\textbf{R}_{2}^{0}=\textbf{0}$
$\displaystyle\textbf{R}_{2}^{k}=\textbf{R}_{1}+\textbf{R}_{1}\textbf{R}_{2}^{k-1}$

where unit ground clauses $r_{1}$ and a recursive clause $r_{2}$ can be
represented as matrices R1, R2 such that $(\textbf{R}_{k})_{i,j}$ = 1 if
$\mathcal{P}\models r_{k}(c_{i},c_{j})$ for constant symbols $c_{i},c_{j}$ and
$(\textbf{R}_{k})_{i,j}$ = 0 otherwise. Addition “+” between two boolean
matrices is defined as (A \+ B)i,j = Ai,j $\lor$ Bi,j. Multiplication
“$\times$” is defined as $(\textbf{A}\times\textbf{B})_{i,j}$ =
$\bigvee_{k=1}^{n}$ Ai,k $\land$ Bk,j and is abbreviated as (AB).

###  3.2 Elementary nets

Figure 2: A one-bounded elementary net (OEN). Places are ellipses. Transitions
are rectangles. Arcs are marked with weights. A black dot is a token. The
initial marking is $m(c_{1})=m(c_{2})=1$. Firing transition $t_{1}$ gives a
new marking $\\{c_{3},c_{4}\\}$.

Petri nets have wide applications in Computer Science
reisig_understanding_2013 . Elementary nets rozenberg_elementary_1998 are a
class of Petri nets. We summarise key concepts and notions of Petri net and
elementary nets from reisig_understanding_2013 ; rozenberg_elementary_1998 . A
Petri net contains two disjoint sets of nodes: places $P$ and transitions $T$.
The flow relation $F\subseteq(P\times T)\cup(T\times P)$ describes arcs
between nodes. All places can be marked by tokens which symbolise the
resources distributed in the network. The assignment of tokens to places is
called a marking (of tokens). For a transition $t$, the pre-set
${}^{\bullet}t$ = $\\{p\in P\,|\,pFt\\}$ contains places whose tokens would be
consumed by $t$ and its post-set $t^{\bullet}$ = $\\{p\in P\,|\,tFp\\}$
contains places that would receive tokens from $t$. Markings represent the
states of a Petri net. The weight of an arc $(p,t)$ or $(t,p)$ is
conventionally written as $\overline{pt}$ or $\overline{tp}$. Elementary nets
are basic Petri nets where transitions consume one token from the pre-set and
relocate a token to the post-set. In an elementary net, an arc weight
$\overline{pt}=1\text{ if }(p,t)\in F$ and $\overline{tp}=1\text{ if }(t,p)\in
F$. An elementary net is a tuple $(P,T,F)$.

A one-bounded elementary net (OEN) is a class of elementary nets that has
places marked with at most one token reisig_understanding_2013 . A marking $m$
in OEN is a characteristic function $m:P\to\\{0,1\\}$ defined as:

$\displaystyle m(x)=\left\\{\begin{aligned} &1\qquad\,\,\,\,\,\text{if a token
is assigned to }x\\\ &0\qquad\,\,\,\,\,\text{otherwise}\end{aligned}\right.$

A place $p$ is marked if $m(p)=1$. An initial marking corresponds to the
initial state of an elementary net and its tokens. Figure 2 shows an example
of OEN with an initial marking of tokens. A transition $t$ is enabled by $m$
if and only if for all arcs $(p,t)\in F$, $m(p)=1$. A step
$m\xrightarrow[]{t}m^{\prime}$ is a fired transition resulting in a new
marking according to:

$\displaystyle m^{\prime}(p)=\left\\{\begin{aligned}
&m(p)-1\qquad\,\,\,\,\,\text{if }(p,t)\in F\\\ &m(p)+1\qquad\,\,\,\,\,\text{if
}(t,p)\in F\\\
&m(p)\qquad\qquad\,\,\,\,\,\,\text{otherwise}\end{aligned}\right.$ (1)

A marking $m^{\prime}$ is reachable from a marking $m$ if and only a finite
step sequence exists. A place is reachable from $m$ if and only if
$m^{\prime}$ is a reachable marking from $m$ and $m^{\prime}(p)=1$. A marking
can be written as a row vector ${[}1,0,1,...{]}$ where each element indicates
if the respective place has a token. A transition can also be represented by a
row vector ${[}z_{1},z_{2},...,z_{k}{]}$ where $z_{i}$ reflects the assignment
of tokens at the pre-set and post-set places. The consecutive firings of
transition update markings and the global state of the system. This represents
the dynamic changes in the distribution of tokenised resources.

##  4 Representing a GEM as a datalog program

###  4.1 Modelling a GEM with an OEN

|
---|---

(a)

(b)

Figure 3: OENs representing the qualitative relationship “$H_{2}$ and $O_{2}$
produce $H_{2}O$”. (a) Each place has up to one token. Each token represents
the abundance of a chemical and this marking means $H_{2}$ and $O_{2}$ are
available. (b) This model qualitatively describes that if $H_{2}$ and $O_{2}$
are available, then $H_{2}O$ becomes available by firing the transition.

We discuss the relationship between GEMs and OENs. Petri nets are commonly
used formalism for representing complex biological systems such as genome-
scale metabolic network models and diagnosing their dynamic behaviours
sahu_advances_2021 . A genome-scale metabolic network model contains
biochemical reactions of substances. Reactants are substances initially
involved in a reaction and the products are the chemicals yielded by a
reaction. A reaction is reversible when the conversion of reactants to
products and products to reactants can take place simultaneously. Reversible
reactions are commonly treated as two reactions in opposite directions. A
reaction involves chemical substances $x_{i}$ and $y_{j}$:

$\displaystyle\text{Irreversible: }x_{1}+x_{2}+...+x_{m}$
$\displaystyle\longrightarrow y_{1}+y_{2}+...+y_{n}$
$\displaystyle\text{Reversible: }x_{1}+x_{2}+...+x_{m}$
$\displaystyle\longleftrightarrow y_{1}+y_{2}+...+y_{n}$

A GEM is an OEN in the sense that a reaction is a transition, and reactants
and products are places. Figure 3 presents a basic OEN modelling the water
synthesis reaction involving $H_{2}$, $O_{2}$ and $H_{2}O$.

Given an auxotrophic mutant experiment, we model a GEM as an OEN to find the
set of synthesisable substances by following the metabolic reaction pathways.
This requires computing the necessary set of places for transitions to fire.
Transitions need to be enumerated to examine if their necessary conditions can
ever be met for them to fire. Every simulation of GEM involves a definition of
the growth medium consisting of essential chemicals for cellular metabolism
such as amino acids, inorganic salts and carbon sources. The definition of a
growth medium serves as the basis of metabolic processes and gives an initial
marking of OEN. Reachable places extending this marking are additional
synthesisable substances.

###  4.2 Transforming an OEN into datalog

Each auxotrophic mutant experiment is treated as a datalog query to find all
reachable places based on an initial marking. We evaluate queries of OEN
reachability by transforming the OEN into a linear and immediately recursive
datalog program. A datalog program consisting of clauses with arity two can be
considered a directed graph where edges joining nodes represent the relation
between terms. However, the multiplicity of the involved entities in a
relation cannot be expressed this way. Instead, places and transitions are
represented as constant and predicate symbols. Each argument of our
transformed datalog program is a concatenation of constant symbols, e.g.
$h2\\_o2$. For our specific purpose, since all place nodes in the OEN are
related to metabolites, every new concatenated constant symbol describes a set
of metabolic substances.

This transformation with constant symbolic concatenation generalises from a
directed graph into a hypergraph berge_hypergraphs_1989 , a directed graph
where an edge can connect to any number of nodes. Formally, a hypergraph
$G=(V,E)$ is a tuple containing a set of nodes $V$ and a set of edges
$E\subseteq 2^{V}\times 2^{V}$ between arbitrary subsets of $V$. Hypernodes
are sets of nodes connected by such edges. The incidence graph of a hypergraph
is bipartite in that for $v\in V$ and $e\in E$, $e$ is connected to $v$ if and
only if $v$ and $e$ are nodes connected by an edge in the bipartite graph.
Given an OEN, we consider a hypergraph with an incidence graph of the same
bipartite structure. A hypernode represents a set of place nodes and is
created for every transition’s input places from the OEN. Hypernodes are
denoted by the new constant symbols that concatenate the member places’ names,
e.g. a hypernode $\\{H_{2},O_{2}\\}$ is represented by the constant symbol
$h2\\_o2$. The number of hypernodes we create is fixed and limited by the
number of transitions. This is because each reaction has a finite number of
reactants and products. The dynamic behaviour of an OEN is specified via an
initial marking. This mapping between an OEN and a hypergraph allows us to
write a datalog program. A recursive clause is added which describes a path in
the hypergraph.

######  Definition 1 (OEN transformation)

The transformation takes an OEN and an initial marking $m$. Let $t$ denote any
transition. Input places of transitions are represented by hypernodes $p^{*}$.
A datalog program $\mathcal{P}$ is written with predicate symbols $t$, $r_{1}$
and $r_{2}$ and variables X, Y and Z:

  * •

write $token(m,p^{*})$ if the marking $m$ includes all members of a hypernode
$p^{*}$

  * •

write $token(m,q)$ if the marking $m$ includes a place $q$

  * •

write $t(p_{1}^{*},p_{2}^{*})$ if $t$ connects two hypernodes $p_{1}^{*}$ and
$p_{2}^{*}$

  * •

write $t(p^{*},q)$ if $t$ connects a hypernode $p^{*}$ and a place $q$

  * •

write $r_{1}(X,Y)\leftarrow token(X,Y)$ and $r_{1}(X,Y)\leftarrow t(X,Y)$ for
every $t$ clause written

  * •

write $r_{2}(X,Y)\leftarrow r_{1}(X,Y)$ and $r_{2}(X,Y)\leftarrow
r_{1}(X,Z),r_{2}(Z,Y)$

The transformed program is linear and immediately recursive since it contains
only one recursive clause and the recursive predicate appears only once in the
head of this clause. We illustrate how reachable places from an initial
marking can be computed by evaluating the transformed program.

###### Example 1

Let the OEN in Figure 2 represent a metabolic pathway. Transitions $t_{1}$,
$t_{2}$ and $t_{3}$ are described by $reaction\\_1$, $reaction\\_2$ and
$reaction\\_3$ clauses. We create two hypernodes denoted by $c1\\_c2$ and
$c3\\_c4$. These transitions connect places $\\{c_{3}$, $c_{4}$, $c_{5}\\}$
which are represented by constants $\\{c3$, $c4$, $c5\\}$. For an initial
marking $m_{1}=\\{c_{1},c_{2}\\}$, we can construct a datalog program
$\mathcal{P}_{1}$:
$\displaystyle\mathcal{P}_{1}:\left\\{\begin{aligned} token(m1,c1\\_c2).&\\\
token(m1,c1).&\\\ token(m1,c2).&\\\ reaction\\_1(c1\\_c2,c3\\_c4).&\\\
reaction\\_1(c1\\_c2,c3).&\\\ reaction\\_1(c1\\_c2,c4).&\\\
reaction\\_2(c3\\_c4,c5).&\\\ reaction\\_3(c5,c4).&\\\ reaction(X,Y)\leftarrow
token(X,Y).&\\\ reaction(X,Y)\leftarrow reaction\\_1(X,Y).&\\\
reaction(X,Y)\leftarrow reaction\\_2(X,Y).&\\\ reaction(X,Y)\leftarrow
reaction\\_3(X,Y).&\\\ pathway(X,Y)\leftarrow reaction(X,Y).&\\\
pathway(X,Y)\leftarrow reaction(X,Z),&\\\
pathway(Z,Y).&\end{aligned}\right\\}$

The right-hand-side figures are the OEN from from Figure 2 and the hypergraph
represented by $\mathcal{P}_{1}$. Both graphs show that the place $c_{2}$ is
immediately reachable from the initial marking $m_{1}=\\{c_{1},c_{2}\\}$. The
hypernode $\\{c_{3},c_{4}\\}$ is reachable from via the transition $t_{1}$.
The place $c_{5}$ is reachable in two steps via transitions $t_{1}$ and
$t_{2}$. We can compute these by evaluating groundings of $pathway(m1,Y)$ in
the Herbrand model $M(\mathcal{P}_{1})$ which gives us $pathway(m1,c2)$,
$pathway(m1,c3\\_c4)$ and $pathway(m1,c5)$.

We prove the correctness of computing reachable places in an OEN by evaluating
the transformed datalog program. We first show the effect of a step (a fired
transition) on updating markings. A step in an OEN updates the set of marked
places.

###### Lemma 1

Let $OEN=(P,T,F)$. For a transition $t\in T$, marking $m$ and $m^{\prime}$,
$\land_{p_{i}\in^{\bullet}t}m(p_{i})=1\to\land_{p_{j}\in
t^{\bullet}}m^{\prime}(p_{j})=1$.

###### Proof

It follows from Equation 1 that all tokens at the places in the pre-set of a
transition are consumed and places in the post-set of a transition are
assigned tokens.

Lemma 1 shows that when a transition $t$ fires, in the new marking, tokens are
re-assigned to other places. Given a marking, the set of marked places is a
subset of $P$. The set of all token-to-place assignments is a power set
$\Omega=2^{P}$ and it is the union of all possible markings. From the property
shown in Lemma 1, we define an operator called step consequence
$\mathcal{T}_{F}:\Omega\to\Omega$.

###### Definition 2

(Step consequence operator $\mathcal{T}_{F}$) Given an OEN = (P, T, F),
$\mathcal{T}_{F}$ is an operator over the power set $\Omega$:

$\displaystyle\forall S\in\Omega,\mathcal{T}_{F}(S)=\\{p_{j}\in P\,|\,p_{j}\in
S^{\prime}\leftarrow\land_{p_{i}\in{{}^{\bullet}}t}p_{i}\in S$
$\displaystyle\text{ where}\,\\{p_{1},...,p_{n}\\}\subseteq S\text{ and
}p_{j}\in t^{\bullet},S^{\prime}\in\Omega\\}$

Definition 3 says that the result of applying $\mathcal{T}_{F}$ to a marking
is a set of places reachable from it. $\mathcal{T}_{F}(\emptyset)$ only
includes the initial marking due to the antecedent of the implication being
false. To obtain the union of all places reachable from an empty marking, we
can apply $\mathcal{T}_{F}$ repeatedly to its results. The recursive
application of $\mathcal{T}_{F}$ can be defined as
$\mathcal{T}_{F}{\uparrow^{0}}=\emptyset$ and
$\mathcal{T}_{F}{\uparrow^{n}}=\mathcal{T}_{F}(\mathcal{T}_{F}{\uparrow^{n-1}})$.
This corresponds to a lattice.

###### Proposition 1

$(\mathcal{T}_{F},\subseteq)$ is a complete lattice.

###### Proof

The proof follows from Definition 2. $(\mathcal{T}_{F},\subseteq)$ is a
complete lattice since $\Omega$ is a set partially ordered by set inclusion
$\subseteq$ where $\mathcal{T}_{F}$ is an operator over the power set
$\Omega$. The top element $P$ is the union of elements in $\Omega$ and the
bottom element $\emptyset$ is the intersection of elements in $\Omega$.

Let $(P,T,F)$ be an arbitrary OEN. Using Proposition 1, we show
$\mathcal{T}_{F}$ over the lattice $(\mathcal{T}_{F},\subseteq)$ has a
fixpoint. We refer to Tarski’s fixpoint theorem tarski_fixpoint , which states
a monotonic operator on a complete lattice has a least fixpoint. We prove that
$\mathcal{T}_{F}$ is a monotonic operator.

###### Lemma 2

$\mathcal{T}_{F}$ is a monotonic operator over $\Omega$. For all
$S_{1},S_{2}\in\Omega$, if $S_{1}\subseteq S_{2}$ then
$\mathcal{T}_{F}(S_{1})\subseteq\mathcal{T}_{F}(S_{2})$.

###### Proof

This follows from Definition 2. We will always extend a marking by adding
places reachable from it. Places in the current marking are reachable and are
retained after applying $\mathcal{T}_{F}$. Therefore, $\mathcal{T}_{F}$ is a
monotonic operator.

We show that the least fixpoint of $\mathcal{T}_{F}$ is the union of reachable
places by proving that $\mathcal{T}_{F}$ is a compact operator.

###### Lemma 3

$\mathcal{T}_{F}$ is a compact operator over $\Omega$. For all
$S_{1}\in\Omega$, if $p\in\mathcal{T}_{F}(S_{1})$ then
$p\in\mathcal{T}_{F}(S_{2})$ for a finite $S_{2}\subseteq S_{1}$.

###### Proof

Assume $p\in\mathcal{T}_{F}(S_{1})$. From Definition 2 and the assumption,
since $\mathcal{T}_{F}$ performs a step, a transition $t$ exists and $S_{1}$
must contain all pre-set places. Take places $S_{2}$ that contains $S_{1}$.
From Lemma 1, we know that firing $t$ give
$S^{\prime}\subseteq\mathcal{T}_{F}(S_{2})$ and $p\in S^{\prime}$.

###### Proposition 2

The operator $\mathcal{T}_{F}$ has a least fixpoint equal to
$\mathcal{T}_{F}{\uparrow^{\alpha}}=\cup_{n\geq
0}\mathcal{T}_{F}{\uparrow^{n}}$.

###### Proof

$\mathcal{T}_{F}$ is a monotonic function. According to Tarski’s fixpoint
theorem tarski_fixpoint , the function $\mathcal{T}_{F}$ has a least fixpoint
$\mathcal{T}_{F}{\uparrow^{\alpha}}$ for some $\alpha\geq 0$. Since
$\mathcal{T}_{F}$ is also a compact operator, if
$p\in\mathcal{T}_{F}(\mathcal{T}_{F}{\uparrow^{\alpha}})$ then
$p\in\mathcal{T}_{F}(\mathcal{T}_{F}{\uparrow^{n}})$ for some $n\geq 0$ and
$\mathcal{T}_{F}{\uparrow^{n}}\subseteq\mathcal{T}_{F}{\uparrow^{\alpha}}$.
The union of subsets of $\mathcal{T}_{F}{\uparrow^{\alpha}}$ is $\cup_{n\geq
0}\mathcal{T}_{F}{\uparrow^{n}}$, which covers every member of the least
fixpoint. Therefore, $\mathcal{T}_{F}{\uparrow^{\alpha}}=\cup_{n\geq
0}\mathcal{T}_{F}{\uparrow^{n}}$.

Proposition 2 shows that when we repeatedly apply the step consequence
operator to the empty set, the least fixpoint is the union of reachable
places. This allows us to compute this least fixpoint by iteratively computing
steps and expanding the set of reachable places. We then show the correctness
of finding reachable places in an OEN via evaluating the transformed linear
and immediately recursive program

###### Theorem 4.1

Let $\mathcal{P}$ be a datalog program transformed from an OEN. For constant
symbols $c_{1}$, $c_{2}$ in $\mathcal{P}$, $\mathcal{P}\models r(c_{1},c_{2})$
if and only if $c_{2}$ is reachable from $c_{1}$ in the OEN given an initial
marking $m$.

###### Proof

Prove the forward implication. Assume $\mathcal{P}\models r(c_{1},c_{2})$.
Constant symbols $c_{1}$ and $c_{2}$ represent place nodes in the OEN. We know
$r(c_{1},c_{2})$ is in the least model. Based on the definition of the
Immediate Consequence Operator van_emden_semantics_1976 ,
$r(c_{1},c_{2})\leftarrow\,B_{1},...,B_{i}$ ($i\geq 0$) is a clause in
$\mathcal{P}$ and
$\\{B_{1},...,B_{i}\\}\subseteq\mathcal{T_{P}}{\uparrow^{\omega}}$. From
Definition 1, we know that $B_{1},...,B_{i}$ are unit ground clauses
describing the pre-sets and post-sets of transitions. There exists $n\geq 0$
such that $c_{1}\in\mathcal{T}_{F}{\uparrow^{n}}$. Following the transitions
$B_{1},...,B_{i}$, we can apply the operator $\mathcal{T}_{F}$ to
$\mathcal{T}_{F}{\uparrow^{n}}$. From Definition 2 and proposition 2, we know
$c_{2}\in\mathcal{T}_{F}{\uparrow^{\alpha}}$ so a set of marked places
$S\in\Omega$ must contain $c_{2}$. Therefore, $c_{2}$ is reachable from
$c_{1}$ given an initial marking $m$.

Now we prove the backward implication. Assume $c_{2}$ is reachable from
$c_{1}$ in the OEN given an initial marking $m$. We know from Proposition 2
that $c_{1},c_{2}\in\mathcal{T}_{F}{\uparrow^{\alpha}}$. From Definition 2,
there exists a finite sequence of steps
$m\xrightarrow[]{t_{1}}m_{1}\xrightarrow[]{t_{2}}...\xrightarrow[]{t_{k}}m^{\prime}$
where $m(c_{1})=m^{\prime}(c_{2})=1$ and $k\geq 0$. From Definition 1, the
transformed program $\mathcal{P}$ contains these transitions as unit ground
clauses $\\{B_{1},...,B_{k}\\}\subseteq\mathcal{B_{P}}$ and recursive clauses
$r(X,Y)\leftarrow r_{1}(X,Y)$ and $r(X,Y)\leftarrow r_{1}(X,Z),r(Z,Y)$. There
is a grounding of the recursive clauses such that $r(c_{1},c_{2})\leftarrow
B_{1},...,B_{k}$. By applying the $\mathcal{T_{P}}$ operator, we know
$\mathcal{P}\models r(c_{1},c_{2})$.

Theorem 4.1 says that a program transformed from an OEN is logically
equivalent to the OEN based on the fixpoint semantics in Proposition 2 for
evaluating OEN reachability queries. In the next section, we explain how we
leverage boolean matrix operations to evaluate datalog transformed from OENs.

##  5 Boolean matrix logic programming

###  5.1 Problem setting

We propose the Boolean Matrix Logic Programming (BMLP) problem. In contrast to
traditional logic program evaluation, BMLP uses boolean matrices to evaluate
datalog programs.

######  Definition 3 (Boolean Matrix Logic Programming (BMLP) problem)

Let $\mathcal{P}$ be a datalog program containing a set of clauses with
predicate symbol $r$. The goal of Boolean Matrix Logic Programming (BMLP) is
to find a boolean matrix R encoded by a datalog program such that
$(\textbf{R})_{i,j}$ = 1 if $\mathcal{P}\models r(c_{i},c_{j})$ for constant
symbols $c_{i},c_{j}$ and $(\textbf{R})_{i,j}$ = 0 otherwise.

We focus on linear and immediately recursive datalog programs with at most two
body literals to simulate a GEM. Our goal is to implement BMLP algorithms that
can compute boolean matrices for evaluating large logic programs representing
GEMs.

####  5.1.1 Creating boolean matrices

From a datalog program $\mathcal{P}$, we define mappings between constant
symbols and binary codes. Consider a ground term $u$ that concatenates $n$
constant symbols $c_{1},c_{2},...,c_{n}$. We can represent $u$ using a $n$-bit
binary code. When $c_{k}$ is in $u$ the $k$-th bit of the binary code is 1,
otherwise, it is 0. For a linear and immediately recursive clause $r$ that
takes such terms as arguments, two boolean matrices $\textbf{R}_{1}$ and
$\textbf{R}_{2}$ are defined by mapping unit ground clauses used by $r$ to
binary codes. The two matrices in logic program format have
$\\{v_{1}(1,b_{1}^{1})$, $v_{1}(2,b_{1}^{2}),...\\}$ and
$\\{v_{2}(1,b_{2}^{1})$, $v_{2}(2,b_{2}^{2}),...\\}$ as their rows where
$v_{1}(i,b_{1}^{i})$ and $v_{2}(i,b_{2}^{i})$ are the $i$-th rows of the
boolean matrices $\textbf{R}_{1}$ and $\textbf{R}_{2}$. In $v_{1}$ and $v_{2}$
clauses, $b_{1}^{i}$ and $b_{2}^{i}$ are binary codes to represent constant
symbols $u_{1}$ and $u_{2}$ as $r$’s first and second arguments.
$\textbf{R}_{1}$ and $\textbf{R}_{2}$ are defined such that
$v_{1}(i,b_{1}^{i})$ and $v_{2}(i,b_{2}^{i})$ hold only when
$\mathcal{P}\models r(u_{1},u_{2})$.

###### Example 2

We show boolean matrices created from $\mathcal{P}_{1}$ in Example 1.
Individual substances are represented by constant symbols $c1$, $c2$ and $c3$.
Sets of chemical substances are denoted by additional concatenated constant
symbols, $c1\\_c2$, $c3\\_c4$. On the right-hand side, we create two boolean
matrices $\textbf{R}_{1}$ and $\textbf{R}_{2}$ from the reaction clauses in
$\mathcal{P}_{1}$. Since the first reaction_1 clause also contains the
constant symbols of the other reaction_1 clauses, to avoid redundancy, only
one binary code is created for these clauses in $\textbf{R}_{1}$ and
$\textbf{R}_{2}$.
$\displaystyle\mathcal{P}_{1}:\left\\{\begin{aligned} ...&\\\
reaction\\_1(c1\\_c2,c3\\_c4).&\\\ reaction\\_1(c1\\_c2,c3).&\\\
reaction\\_1(c1\\_c2,c4).&\\\ reaction\\_2(c3\\_c4,c5).&\\\
reaction\\_3(c5,c4).&\\\ ...&\\\ reaction(X,Y)\leftarrow
reaction\\_1(X,Y).&\\\ reaction(X,Y)\leftarrow reaction\\_2(X,Y).&\\\
reaction(X,Y)\leftarrow reaction\\_3(X,Y).&\\\ pathway(X,Y)\leftarrow
reaction(X,Y).&\\\ pathway(X,Y)\leftarrow reaction(X,Z),&\\\
pathway(Z,Y).&\end{aligned}\right\\}$ $\displaystyle\textbf{R}_{1}\text{
(matrix)}:\begin{bmatrix}0\quad 0\quad 0\quad 1\quad 1\\\ 0\quad 1\quad 1\quad
0\quad 0\\\ 1\quad 0\quad 0\quad 0\quad 0\end{bmatrix}$
$\displaystyle\textbf{R}_{2}\text{ (matrix)}:\begin{bmatrix}0\quad 1\quad
1\quad 0\quad 0\\\ 1\quad 0\quad 0\quad 0\quad 0\\\ 0\quad 1\quad 0\quad
0\quad 0\end{bmatrix}$ $\displaystyle\textbf{R}_{1}\text{
(program)}:\left\\{\begin{aligned} &v_{1}(1,00011).\\\ &v_{1}(2,01100).\\\
&v_{1}(3,10000).\end{aligned}\right\\}$ $\displaystyle\textbf{R}_{2}\text{
(program)}:\left\\{\begin{aligned} &v_{2}(1,01100).\\\ &v_{2}(2,10000).\\\
&v_{2}(3,01000).\end{aligned}\right\\}$

Constant symbols are represented by rows in $\textbf{R}_{1}$ and
$\textbf{R}_{2}$. Constant symbols $c1\\_c2$, $c3\\_c4$ and $c5$ that appear
on the left-hand side of reactions are denoted by binary bits in
$\textbf{R}_{1}$. Constant symbols $c3$, $c4$, $c5$ and $c3\\_c4$ that appear
on the right-hand side of reactions are denoted by binary bits in
$\textbf{R}_{2}$. Each row in the boolean matrices is written as a unit ground
clause. For example, since $\mathcal{P}_{1}\models reaction(c1\\_c2,c3\\_c4)$,
the first row of $\textbf{R}_{1}$ is written a binary code “00011” to
represent $c1\\_c2$ and the first row of $\textbf{R}_{2}$ is a binary code
“01100” to represent $c3\\_c4$. Each “1” binary bit indicates the presence of
a constant symbol.

Each transformed program contains clauses describing an initial marking. These
clauses are represented by a boolean vector v which is a single-row boolean
matrix where the binary bits correspond to the constant symbols.

###### Example 3

We show the mapping of a single row matrix from the token definition from
$\mathcal{P}_{1}$ in Example 1. Individual substances are represented by
constant symbols $c1$ and $c2$. The set of chemical substances is denoted by
an additional concatenated constant symbol, $c1\\_c2$.
$\displaystyle\mathcal{P}_{1}:\left\\{\begin{aligned} token(m1,c1\\_c2).&\\\
token(m1,c1).&\\\ token(m1,c2).&\\\ ...&\\\ reaction(X,Y)\leftarrow
token(X,Y).&\\\ ...&\end{aligned}\right\\}$ $\displaystyle\textbf{v}\text{
(vector)}:\begin{bmatrix}0\quad 0\quad 0\quad 1\quad 1\end{bmatrix}$
$\displaystyle\textbf{v}\text{ (program)}:\left\\{\begin{aligned}
v(1,00011).\end{aligned}\right\\}$

On the right-hand side, the constant symbol $c1\\_c2$ is represented by the
vector v. Since the first token clause contains constant symbols in the other
token clauses, v only includes the right-hand side of the first token clause
to avoid redundancy. In program form, v is written as a unit ground clause
containing a binary code “00011” to represent $c1\\_c2$.

The separation of this vector from $\textbf{R}_{1}$ and $\textbf{R}_{2}$ is
meant to keep $\textbf{R}_{1}$ and $\textbf{R}_{2}$ static. Since the
structure of the OEN does not vary often, we reduce the number of times that
we need to re-create these two matrices while the definition of an initial
marking usually changes across different simulations. We use all three created
matrices, $\textbf{R}_{1}$, $\textbf{R}_{2}$ and v for bottom-up evaluation.

###  5.2 Iterative extension (BMLP-IE)

Algorithm 1 Iterative extension (BMLP-IE)

Input: A $1\times n$ vector v, a $1\times k$ vector t that represents viable
transitions, two $k\times n$ boolean matrices $\textbf{R}_{1}$,
$\textbf{R}_{2}$ that encode unit ground clauses.
Output: Transitive closure $\textbf{v}^{*}.$

1: Let $\textbf{v}^{*}$ = v, $1\times k$ vector $\textbf{v}^{\prime}$ = 0.

2: while $True$ do

3: for $1\leq i\leq k$ do

4: $(\textbf{v}^{\prime})_{i}=(\textbf{t})_{i}$ if (v AND
$(\textbf{R}_{1})_{i,*}$) == $(\textbf{R}_{1})_{i,*}$.

5: end for

6: $\textbf{v}^{*}$ = ($\textbf{v}^{\prime}$ $\textbf{R}_{2}$) + v.

7: $\textbf{v}=\textbf{v}^{*}$ if $\textbf{v}^{*}\neq\textbf{v}$ else break.

8: end while

Given a program $\mathcal{P}$ that contains a linear and immediately recursive
clause $r$, we implement a BMLP algorithm called iterative extension to
perform the bottom-up evaluation of $r$. The iterative extension algorithm
(BMLP-IE) operates on two matrices to iteratively expand terms reachable in a
partially grounded query $r(u,Y)$. Given that $\mathcal{P}$ contains $n$
constant symbols, we represent the constant symbols in $u$ as a $1\times n$
row vector or a $n$-bit binary code $\mathbf{v}$. Based on Section 5.1.1, we
write $k\times n$ boolean matrices $\textbf{R}_{1}$ and $\textbf{R}_{2}$ for
$k$ unit ground clauses used by $r$. $\textbf{R}_{1}$ and $\textbf{R}_{2}$
describe which constant symbols appear first and second in the unit ground
clauses. In addition, a $1\times k$ boolean vector t is used to express
additional control on the viability of transitions in the OEN. Transition can
be turned “on” or “off” via this vector. “On” transitions have the binary
value 1 and “off” transitions have the value 0. This allows us to temporarily
modify the network connectivity without re-creating $\textbf{R}_{1}$ and
$\textbf{R}_{2}$.

We implement BMLP-IE (Algorithm 1) in SWI-Prolog (version 9.2) to compute the
transitive closure of v. In SWI-Prolog, each binary code can be written as an
integer number, e.g. “00011” is the integer “3”. We use the bitwise AND to
locate rows in $\textbf{R}_{1}$ that contain all the 1 elements in v. This
computes indices of unit ground clauses that can extend the set of constant
symbols in v. These constants are added to v via boolean matrix multiplication
and addition. We iterate this process until we find the elements of v no
longer change.

###### Theorem 5.1

Given a $1\times n$ vector v, a $1\times k$ vector t, two $k\times n$ boolean
matrices $\textbf{R}_{1}$ and $\textbf{R}_{2}$, Algorithm 1 has a time
complexity $O(kn^{2})$ for computing the transitive closure $\textbf{v}^{*}$.

###### Proof

In Algorithm 1, the complexity of an iteration in the “while” loop is
$O(k\times n)$ bitwise operations due to the multiplication between a vector
and a $k\times n$ matrix. Each “for” loop pass also performs $O(k\times n)$
bitwise comparisons. Until we find the transitive closure, at least one unit
ground clause needs to be added in each “while” loop iteration. Therefore,
there are at most $n$ iterations which require $O(k\times n^{2})$ bitwise
operations. $\square$

Theorem 5.1 shows that BMLP-IE has a polynomial runtime cubic in the size of
the boolean matrices. Evaluating the program $\mathcal{P}$ using SLD
resolutions has the same time complexity and the runtime of BMLP-IE primarily
benefits from bit-wise operations. For SLD resolutions, the size of the subset
of the Herbrand base provable from a query $r(u,Y)$ is bounded by $O(n)$. Each
member of this subset can be found via a depth-first search. Since there are
$k$ unit ground clauses and each variable can be grounded by at most $n$
constants, this query requires at most $n$ search in a tree of $O(k\times n)$
nodes. For OEN reachability queries, it is important to note that the size of
boolean matrices grows linearly with the number of hypernodes because they are
both bounded by the number of transitions.

##  6 Active learning

###  6.1 Hypothesis compression

We use the compression score of a hypothesis $h$ to optimise the posterior
probability. A learner can derive the relative frequency $p(h)$ for $h$ to be
chosen based on the encoding of $h$. In an optimal binary bit coding scheme
shannon_mathematical_1963 , the length of $h$ denoted by $size(h)$ is
$-log_{2}(p(h))$. The compactness of $h$’s encoding reflects its prior
probability $p(h)=2^{-size(h)}$. We can use the compression function to
compute the prior $p(h)$. Given a set of training examples $E$, the
compression of a hypothesis progol ; ase_progol can be interpreted as

$\displaystyle compression(h,E)=size(E)-MDL_{h,E}$

In addition, the compression relates to the posterior probability of $h$
progol :

$\displaystyle\frac{p(h|E)}{p(E|E)}=2^{compression(h,E)}$

We can re-arrange this to obtain the normalised posterior probability
ase_progol :

$\displaystyle p^{\prime}(h|E)=\frac{2^{compression(h,E)}}{\sum_{h_{i}\in
H}2^{compression(h_{i},E)}}$

The posterior probability is maximal when compression is maximal. By searching
for a hypothesis with maximum posterior probability, a learner can maximise
its expected predictive accuracy progol . To learn an accuracy hypothesis, the
goal is to find a hypothesis with the highest compression score.

The compression function in $BMLP_{active}$ follows the Minimal Description
Length (MDL) principle conklin_complexity-based_1994 . The most probable
hypothesis $h$ for the training examples $E$ should be compact and have the
fewest disagreements between $E$ and predictions made by $h$. This minimises
$L(h)+L(E|h)$, where $L(h)$ is the descriptive complexity of $h$ and $L(E|h)$
is the descriptive complexity of encoding $E$ using $h$. We consider a
compression function from ase_progol based on the MDL principle:

$\displaystyle
compression(h,E)=|E^{+}|-\frac{|E^{+}|}{pc_{h}}(size(h)+fp_{h})$

where $E^{+}$ is the set of positive examples seen by the active learner,
$pc_{h}$ is the positive coverage by $h$ and $fp_{h}$ is false positives
covered by $h$. In other words, this compression function favours a general
hypothesis with high coverage and penalises over-general hypotheses that
incorrectly predict negative examples. The generality of $h$ can be estimated
by computing its coverage $pc_{h}$ for a set of unlabelled instances progol .
To determine the hypotheses’ compression scores, $BMLP_{active}$ looks at all
available experiment instances to obtain a good estimate of a hypothesis’s
generality. Learning with increasingly large hypothesis space and instance
space thus requires efficient inferences in the background knowledge. Our BMLP
approach addresses this by manipulating boolean matrices for bottom-up
evaluations.

###  6.2 Expected cost of experiments

Experiment planning is often restricted by time and other resources. A rule of
thumb for experiment design is to achieve optimal empirical outcomes within a
finite budget. The problem of designing experiments to support or refute
scientific hypotheses involves a series of binary decisions. This process is
analogous to playing the game of “Guess Who” where players ask a sequence of
“yes” or “no” questions whose answers binarily partition the candidate
hypotheses. The outcome of an experiment $t$ splits hypotheses $H$ into
consistent hypotheses $H_{t}$ and inconsistent hypotheses $\overline{H_{t}}$.
The path of selection and outcomes of experiments can be represented as a
binary decision tree. To construct an optimal tree, an active learner should
seek to label instances consistent with up to half of the hypotheses
mitchell_generalization_1982 .

In addition to this principle, we consider a user-defined experiment cost
function $C_{t}$ of experiment $t$. The experiment selection is represented by
a binary tree with paths annotated by the costs of experiments. The optimal
selection is therefore a binary tree with the minimum expected overall cost.
The minimum expected cost of performing a set of candidate experiments $T$ can
be recurrently defined as:

$\displaystyle EC(\emptyset,T)=0$ $\displaystyle EC({h},T)=0$ $\displaystyle
EC(H,T)=min_{t\in
T}\,[C_{t}+p(t)EC(H_{t},T-{t})+(1-p(t))EC(\overline{H_{t}},T-{t})]$

To estimate this recurrent cost function, $BMLP_{active}$ uses the following
heuristic function from ase_progol to approximate optimal cost selections:

$\displaystyle EC(H,T)\approx min_{t\in T}\,[C_{t}$
$\displaystyle+p(t)(mean_{t^{\prime}\in T-\\{t\\}}C_{t^{\prime}})J_{H_{t}}$
$\displaystyle+(1-p(t))(mean_{t^{\prime}\in
T-\\{t\\}}C_{t^{\prime}})J_{\overline{H_{t}}}]$

where $p(t)$ is the probability that the outcome of the experiment $t$ is
positive and $J_{H}=-\sum_{h\in H}p(h)\,log_{2}(p(h))$. The probability $p(h)$
is the prior of hypothesis $h$ calculated from the compression function. To
estimate $p(t)$, we compute the sum of the probabilities of the hypotheses
consistent with a positive outcome of $t$. While we only focus on auxotrophic
mutant experiments in this work, users can flexibly define the experiment cost
function for different types of experiments.

###  6.3 Active learning sample complexity

For some hypothesis space $H$ and background knowledge $BK$, let $V_{s}$
denote the version space of hypotheses consistent with $s$ training examples.
For an active version space learner that selects one instance per iteration,
$|V_{s}|$ denotes the size of the version space at the iteration $s$ of active
learning. The shrinkage of the hypothesis space can be represented by the
reduction ratio $\frac{|V_{s+1}|}{|V_{s}|}$ after querying the label of the
$s$-th instance hocquette_how_2018 . The minimal reduction ratio
$p(x_{s},V_{s})$ is the minority ratio of the version space $V_{s}$
partitioned by an instance $x_{s}$.

###### Definition 4

(Minimal reduction ratio hocquette_how_2018 ) The minimal reduction ratio over
the version space $V_{s}$ by sampled instance $x_{s}$ is

$\displaystyle p(x_{s},V_{s})=\frac{min(\,|\\{{h\in V_{s}|\,h\cup BK\models
x_{s}\\}|},|\\{h\in V_{s}|\,h\cup BK\not\models x_{s}\\}|)}{|V_{s}|}$

The minimal reduction ratio of an actively selected instance can be computed
before it is labelled. While in reality training examples might not be as
discriminative, the optional selection strategy is to select instances with
minimal reduction ratios as close as possible to $\frac{1}{2}$ with the
ability to eliminate up to 50% of the hypothesis space. We describe the sample
complexity advantage of active learning over random example selection by the
following bound of an active version space learner’s sample complexity. A
passive learner is considered a learner using random example selection since
it does not have control over the training examples it uses.

######  Theorem 6.1 (Active learning sample complexity bound)

For some $\phi\in[0,\frac{1}{2}]$ and small hypothesis error $\epsilon>0$, if
an active version space learner can select instances to label from an instance
space $\mathcal{X}$ with minimal reduction ratios greater than or equal to
$\phi$, the sample complexity $s_{active}$ of the active learner is

$\displaystyle s_{active}\leq\frac{\epsilon}{\epsilon+\phi}s_{passive}+c$ (2)

where c is a constant and $s_{passive}$ is the sample complexity of learning
from randomly labelled instances.

###### Proof

We consider an arbitrary distribution of the minimal reduction ratio over
$\mathcal{X}$ bounded by $[0,\frac{1}{2}]$. Suppose a version space learner
can select instances with a minimal reduction ratio greater than $\phi$. Since
$\phi$ represents the ratio of the minority partition, the number of candidate
hypotheses remaining in the version space $V_{s}$ after $s$ labelled instances
is bounded by the size of the majority partitions $(1-\phi)^{s}|H|$ where $H$
is the hypothesis space. The probability that any hypothesis in $V_{s}$ with
an error larger than $\epsilon$ is consistent with $s$ training examples is:

$\displaystyle p_{error}$ $\displaystyle=(1-\epsilon)^{s}|V_{s}|$
$\displaystyle\leq(1-\epsilon)^{s}(1-\phi)^{s}|H|$

$p_{error}$ should be small so suppose for small $\delta>0$,

$\displaystyle(1-\epsilon)^{s}(1-\phi)^{s}|H|=\delta$ (3)

Equation 3 corresponds to the probability of the extreme case when only the
minority of the remaining hypotheses can be eliminated by each labelled
instance. By rearranging Equation (3), the number of training examples
required by this learner in this extreme case is:

$\displaystyle
s=\frac{1}{-ln(1-\epsilon)-ln(1-\phi)}(ln(|H|)+ln(\frac{1}{\delta}))$

For an arbitrary active learner whose membership queried instances also have a
minimal reduction ratio greater than or equal to $\phi$, the number of
actively selected instances $s_{active}$ is upper bounded by $s$ for reducing
the size of version space to $|V_{s}|$. Since $(1-u)<e^{-u}$ for small $u>0$,
it is true that $-ln(1-u)>u$. So the following holds:

$\displaystyle
s_{active}<\frac{1}{\epsilon+\phi}(ln(|H|)+ln(\frac{1}{\delta}))$

According to the Blumer bound blumer_occams_1987 , for a version space learner
that randomly samples examples to produce a hypothesis of error at most
$\epsilon$ with a probability less than or equal to $\delta$, the sample
complexity is $s_{passive}>\frac{1}{\epsilon}(ln(|H|)+ln(\frac{1}{\delta}))$.
Following this, it must be true that
$\epsilon\,s_{passive}>(\epsilon+\phi)\,s_{active}$. So we can derive the
Inequality 2 in the theorem:

$\displaystyle s_{active}\leq\frac{\epsilon}{\epsilon+\phi}s_{passive}+c$

where c is a constant. $\square$

This theorem shows the sample complexity gain of active learning over learning
from randomly sampled instances. Theorem 6.1 says that the number of instances
needed to select by active learning should be some factor smaller than the
number of randomly sampled examples given some predictive accuracy level.
Given a desirable predictive error $\epsilon$, this factor can be computed if
the minimal reduction ratio $\phi$ can be estimated from the instance space.
When selected instances have a larger minimal reduction ratio, the sample
complexity gain would be more significant. On the other hand, if instances
have low discriminative power and the minimal reduction ratio tends to zero,
the sample complexity gain would be minimal. In Section 8, we empirically show
empirical evidence supporting this theorem for $BMLP_{active}$ and randomly
sampled experiments when learning gene functions in a GEM.

##  7 Implementation

###  7.1 Encoding a metabolic network

Our encoding of a metabolic network is based on the Robot Scientist
King04:RobotScientist ; ase_progol . The Robot Scientist is implemented using
the ILP system Progol progol and finds the synthesisable set of metabolites by
traversing a metabolic network. It uses a metabolic network as a deductive
database and answers a query of user-defined gene deletion and experimental
conditions based on SLD resolutions. Its deductive database is written in
Prolog syntax and uses Prolog list operations to manipulate constant symbols
representing metabolic substances. Metabolic reactions are enumerated based on
the set of substances currently synthesisable. The products of a metabolic
reaction extend the list of synthesisable substances. Since the number of
metabolic reactions is finite, the closure of members in this list needs to
contain essential substances for the growth of an organism. In an auxotrophic
mutant experiment outcome prediction, a gene function is removed. Some key
metabolic reaction pathways might no longer be viable due to the removal,
leading to insufficient substances for growth. In this case, the predicted
experiment’s outcome is a lack of growth and a phenotypic effect.

In the Robot Scientist, the metabolic network is defined by $enzyme/6$
clauses. Each $enzyme/6$ clause contains an artificial identifier for an
enzyme, gene locus identifiers, reaction identifiers, a reaction reversibility
indicator and reactant and product identifiers. A gene locus corresponds to a
physical site within a genome and is uniquely associated with a gene
identifier. An artificial enzyme identifier is created for each unique gene
locus identifier. An auxotrophic mutant experiment is described as a query of
$phenotypic\\_effect/2$ whose outcome is a binary prediction. A
$phenotypic\\_effect/2$ query has a gene identifier and a list of reagents
which are optional chemical substances that provide additional nutrients to
facilitate metabolic reactions. These reagents are added to a list of pre-
defined base growth medium compounds. Gene deletions are realised via
$codes/2$ clauses where an enzyme corresponding to the deleted gene would be
removed from the set of metabolic reactions to be enumerated.

###### Example 4

In the following program $\mathcal{P}_{2}$, the $codes/2$ clause says that the
gene locus “b2379” is responsible for the enzyme “enz201”. The $enzyme/6$
clause defines the reversible reaction “alata_l” (L-alanine transaminase).
This reaction involves reactants “akg_c” (Oxoglutarate), “ala__L_c”
(L-alanine) and products “glu__L_c” (L-glutamate), “pyr_c” (Pyruvate). An
integer value in {1,2} indicates the reaction reversibility and the integer
“2” says that the reaction is reversible.

$\displaystyle\mathcal{P}_{2}:\left\\{\begin{aligned} codes(b2379,enz201).&\\\
enzyme(enz201,[b2379],[alata\\_l],2,[[akg\\_c,ala\\_\\_L\\_c]],&\\\
[[glu\\_\\_L\\_c,pyr\\_c]]).&\end{aligned}\right\\}$

The following $phenotypic\\_effect/2$ query describes an auxotrophic mutant
experiment where the gene “b2379” is removed and an additional metabolite
“glyc_e” (Glycerol) is introduced in the growth medium.

$\displaystyle?-phenotypic\\_effect(b2379,[glyc\\_e]).$

iML1515 separates metabolites into different compartments in the cell. The
suffices “e” (extracellular space) and “c” (cytosol) represent two
compartments, which divide the metabolic network into sub-networks.

Our encoding of the GEM iML1515 also contains $codes/2$ and $enzyme/6$
definitions. These clauses enable us to describe metabolic reactions via
boolean matrices. We convert gene-reaction-metabolite mappings and standard
identifiers in iML1515111The iML1515 model is described in the text-based
JavaScript Object Notation (JSON). It uses standard identifiers from the
Biochemical Genetic and Genomic (BiGG) king_bigg_2016 knowledge base of large-
scale metabolic reconstructions. into a logic program containing $codes/2$ and
$enzyme/6$ clauses. iML1515 uses standard nomenclature so the logic program
encoding can be applied to other genome-scale metabolic networks. The logic
program encoding can also be linked to publications and external databases to
obtain further details on genes, proteins, reactions and metabolites.

###  7.2 Simulating a GEM using BMLP-IE

Algorithm 2 BMLP phenotype prediction

Input: Background knowledge $BK$, an experiment query $q$.
Output: A binary phenotypic effect classification: 1 represents a phenotypic
effect and 0 denotes no phenotypic effect.

1: From $reaction/2$ clauses in $BK$, create boolean matrices
$\textbf{R}_{1}$, $\textbf{R}_{2}$.

2: From biomass metabolites in $BK$, create vector $\textbf{v}_{biomass}$.

3: From growth medium metabolites in $q$, create vector v.

4: From knockout genes in $q$, create vector t for viable reactions.

5: return 1 if (BMLP-IE(v, t, $\textbf{R}_{1}$, $\textbf{R}_{2}$) AND
$\textbf{v}_{biomass})\neq\textbf{v}_{biomass}$ else return 0.

To simulate the GEM iML1515 using BMLP-IE, $reaction/2$ clauses are created
from $enzyme/6$ clauses. Each clause contains reactants and product
identifiers. Each reversible reaction is defined as two clauses where
reactants and products are swapped. One difference between $BMLP_{active}$ and
the Robot Scientist is that $BMLP_{active}$ allows multiple gene deletions,
e.g. $phenotypic\\_effect([b2290,b2379]$, $[glyc\\_e])$. This enables the
learning of digenic interactions such as gene-isoenzyme functions. Since
multiple genes can be responsible for the same enzyme, single knockouts may
not have phenotypic effects. A hypothesis of a gene-isoenzyme function
contains one $codes/2$ clause and one $enzyme/6$ clause. The $enzyme/6$ in the
hypothesis would be a replicate of an existing clause in the background
knowledge with a new enzyme identifier.

Another difference is that boolean matrices are created using the $reaction/2$
clauses for bottom-up evaluation. Algorithm 2 shows the application of BMLP-IE
(Algorithm 1) for finding the transitive closure of synthesisable substances.
The two boolean matrices inputs to the BMLP-IE algorithm represent reactants
and products participating in the reactions. The growth medium condition is
encoded as the input vector to the BMLP-IE algorithm. This vector would be
extended to represent substances that are synthesisable by the organism in an
auxotrophic mutant experiment. Specifying knockout genes excludes a subset of
reactions. A vector computed during an intermediate step in the BMLP-IE
algorithm denotes viable reactions given the set of knockout genes. A viable
reaction is indicated by a binary bit 1 in this vector. Blocked reactions are
excluded based on $codes/2$ and $enzyme/6$ definitions. In BMLP-IE, boolean
operations compute the transition between metabolic states and return the
closure set of synthesisable substances. A binary phenotypic classification is
made according to the inclusion of essential substances in the closure set,
indicating cell growth comparable to an unedited wild-type MG1655 E. coli
strain. A binary classification of “1” corresponds to a phenotypic effect and
“0” is no phenotypic effect. The essential substances come from the biomass
function of the wild-type in iML1515 iML1515 222We use the wild-type biomass
function specified by iML1515 from http://bigg.ucsd.edu/ king_bigg_2016 . We
use the qualitative version of this biomass function, meaning that we only
care about the substances in the biomass, not their reaction coefficients..
The biomass function describes the necessary substance composition for the
growth of the organism.

###  7.3 Active learning using BMLP

Algorithm 3 $BMLP_{active}$

Input: Hypotheses $H$, background knowledge $BK$, experiment instances $T$,
experimental cost budget $C$, total number of experiment $N$ and oracle $O$.
Output: A hypothesis $h_{max}\in H$.

1: Let $T_{selected}=\emptyset$, $V_{0}=H$.

2: while $cost(T_{selected})\leq C$ and $|V_{j}|>1$ and $|T_{selected}|<N$ do

3: if $T_{selected}==\emptyset$ then

4: Select $t_{j}$ with the largest minimum reduction ratio $p(t_{j},V_{j})$

5: from remaining experiments in $T$ to add to $T_{selected}$.

6: else

7: Add experiment $t_{j}\in T$ that has the minimum cost
$EC(V_{j},T^{\prime})$ to $T_{selected}$.

8: end if

9: break if $cost(T_{selected})>C$

10: Consult the label $l$ of $t_{j}$ from the oracle $O$.

11: $\textbf{R}_{i,j}$ = $BMLP\\_phenotype\\_prediction$($h_{i}\cup BK$,
$t_{j}$) for all $h_{i}\in V_{j}$

12: $V_{j+1}=\\{h_{i}\,|\,h_{i}\in V_{j}\text{ and }\textbf{R}_{i,j}=l\\}$.

13: For $h_{i}\in V_{s+1}$, compute $compression(h_{i},T_{selected})$.

14: $h_{max}$ is the hypothesis with the largest compression score.

15: end while

Active learning in $BMLP_{active}$ is described by Algorithm 3. The input
background knowledge is created from a GEM model. The user needs to provide a
set of candidate experiments, candidate hypotheses and an experiment cost
function. Labels of experiment instances are requested from a data source,
e.g. a laboratory or an online dataset. Labelled experimental data are
considered ground truths. Initially, we randomly shuffle all candidate
experiments and select an experiment with the most discriminative power if its
cost is under budget. Alternatively, a discriminative experiment with the
lowest cost can be selected for a fixed initialisation.

In each active learning cycle in $BMLP_{active}$, hypotheses inconsistent with
ground truth experimental outcomes are pruned. $BMLP_{active}$ selects
experiments to minimise the expected value of a user-defined cost function
(see Section 6). We apply BMLP-IE to obtain predicted experiment outcomes and
compare them with experimental data. For a hypothesis space $H$ and a set of
experiment instances $T$, BMLP-IE produces predicted outcomes in a boolean
matrix R of size $|H|\times|T|$. Given the j-th actively selected instance
$t_{j}\in T$ and $h_{i}\in V_{j}\subseteq H$, $(\textbf{R})_{i,j}=1$ if
$h_{i}\cup BK\models t_{j}$ and otherwise $(\textbf{R})_{i,j}=0$. We refer to
R as a classification table. In experiments from Section 8, we construct
classification tables by predicting phenotypic effects for all combinations of
candidate experiments and hypotheses. This avoids repeating simulations in
hypothesis pruning and compression calculation. When computing hypothesis
compression, we estimate the generality of a hypothesis by calculating its
coverage for all unlabelled instances based on the classification table.
Currently, we only consider labels from synthetic or online phenotype data.
However, $BMLP_{active}$ can be coupled with a high-throughput experimental
workflow to automate experiments.

##  8 Experiments

Since a logic program representing a GEM describes many biological entities
and intricate relations between them, abduction of gene function hypotheses
and in-vivo validations requires exploring a large hypothesis space and
instance space. Both are challenging computationally and empirically.
$BMLP_{active}$ aims to address these by evaluating a large logic program with
boolean matrices and actively selecting experiment instances to perform. In
our experiments333All experiments are performed on a workstation with Intel(R)
Core(TM) i9-7900X CPU 3.30GHz (20 CPUs). , we intend to demonstrate that the
$BMLP_{active}$ framework can efficiently perform phenotypic predictions and
actively suggest informative experiment instances that can reduce experiment
costs.

In this work, we define the experimental cost function as the cost of the
reagents which are optional nutrient substances used in selected experiments.
We also assume when added to growth media, reagents are used in the same
volume. In addition, we do not explicitly account for plastic consumables,
equipment deprecation or specify the time to carry out each experiment, which
we consider fixed and are affected mainly by the number of experiments. We
also do not account for costs incurred due to different experiment types such
as different measurement timescales that prevent experiments from being
performed concurrently. To maintain consistent scaling for comparisons, the
total experimental reagent cost is the amount of unit cost required by
selected experiments. We consider the unit cost of the reagents as the cost of
the cheapest reagent. We calculate the total experimental reagent cost as the
sum of the reagent costs444We provide the costs of reagents in the Appendix.
in selected experiments divided by the unit cost.

Q1: Can BMLP-IE reduce the runtime of simulating GEM?

Using only Prolog to simulate metabolic networks requires grounding viable
reaction pathways based on SLD resolutions. We aim to demonstrate that a
general-purpose Prolog interpreter can be augmented with our BMLP approach for
faster phenotypic simulations given logically encoded metabolic network
models. $BMLP_{active}$ extends the logical representation of a metabolic
network in the Robot Scientist King04:RobotScientist ; ase_progol by
additionally applying boolean matrix computation via BMLP. Therefore, we
consider the Robot Scientist’s logical encoding as a benchmark. Since the
Robot Scientist was implemented using Progol progol , we adapt its metabolic
network encoding to SWI-Prolog. To answer Q1, we compare the simulation
runtime of $BMLP_{active}$, which uses SWI-Prolog and BMLP-IE, against the
Robot Scientist’s encoding in SWI-Prolog.

Q2: Does BMLP-IE accurately predict phenotypes?

BMLP-IE makes binary phenotypic predictions and accounts for the necessary
growth conditions of the organism. In contrast to qualitative predictions
produced by BMLP-IE, FBA palsson_systems_2015 is more commonly used in the
literature to quantitatively predict the cell growth rate based on fluxes of
biomass metabolites at the steady state. We compare BMLP-IE phenotypic
predictions against those given by FBA from the iML1515 model in iML1515 . To
answer Q2, the comparison between BMLP-IE and FBA predictions is made with
respect to ground truth experimental data also in iML1515 .

Q3: Can $BMLP_{active}$ reduce the cost for learning gene functions in GEM?

We consider two measures of experiment cost: the experimental resource cost
and the number of experiments. These two costs are related since the number of
experiments generally is the main driver of experimental cost. We aim to show
that experiments selected by $BMLP_{active}$ reduce overall experimental
reagent cost compared to random experiment selection. We additionally
demonstrate that $BMLP_{active}$ lowers the number of experiments (the sample
complexity) needed to learn accurate gene functions. Gene function recovery
experiments are conducted to answer Q3. To demonstrate that $BMLP_{active}$
can recover known gene functions, we delete parts of iML1515 to create an
incomplete knowledge base artificially. Active learning from the incomplete
model is a mock of an actual discovery. Active learning aims to recover
deleted gene functions within a reagent cost budget and with a limited number
of experiments. We also examine Q3 for recovering isoenzyme functions which
require navigating a large experimental design space.

Q4: Does Theorem 6.1 correctly describe the relationship between the sample
complexity of $BMLP_{active}$ and random experiment selection?

Theorem 6.1 shows the advantage of active learning in improving the efficiency
of experimentation. We empirically examine the sample complexity bound
specified in Theorem 6.1. Let the sample complexity of $BMLP_{active}$ and a
random experiment sampling learner be $s_{active}$ and $s_{passive}$
respectively. From Theorem 6.1 we can obtain the following inequality:

$\displaystyle\frac{\epsilon}{s_{active}}\geq\frac{\epsilon+\phi}{s_{passive}}+c\text{
with c as constant}$

In order to answer Q4, we use $BMLP_{active}$ and random experiment sampling
to recover gene functions as described above. We calculate the ratios of the
error rates of learned gene functions over the number of experiments selected.

Q5: Can $BMLP_{active}$ learn digenic interactions?

While isoenzyme functions are identified as one source of errors in iML1515
bernstein_evaluating_2023 , full exploration of the gene-isoenzyme relations
is laborious. This motivates us to answer Q5. We perform gene function
recovery experiments using $BMLP_{active}$ to rediscover a key isoenzyme for
demonstration. We focus on the aromatic amino acid biosynthesis pathway. It is
a critical pathway for micro-organism survival that contains important gene
targets for engineering high aromatic amino acid production strains. Although
this pathway has been intensively studied for decades, gaps remain in the
understanding of gene-gene relations price_filling_2018 . We show that an
isoenzyme function in the aromatic amino acid pathway can be accurately
recovered via active learning. This demonstrates when $BMLP_{active}$ is
integrated with an experimental workflow that systematically perturbs genes
under different growth conditions, we can potentially reveal novel insights
for metabolic engineering and strain optimisation.

###  8.1 Experiment 1: BMLP-IE simulations and predictions

| Single thread (seconds) | 20 CPUs (seconds)
---|---|---
SWI-Prolog | $37.842\pm 9.668$ | $5.089\pm 0.119$
SWI-Prolog + BMLP-IE | $\mathbf{0.220\pm 0.099}$ | $\mathbf{0.061\pm 0.008}$

Table 1: Simulation runtimes in CPU time. BMLP-IE leads to a 170 times
improvement in simulation time efficiency. Multi-threading BMLP-IE enhances
runtime efficiency 600 times compared to base SWI-Prolog.

Materials. In each simulation, we used $BMLP_{active}$ to predict a phenotype
in an auxotrophic single-knockout mutant experiment. Every experiment removes
one gene from the biological system and uses a growth medium. The GEM iML1515
has been validated against 16 carbon sources as optional nutrients in iML1515
. We focused on the same experimental conditions for simulations with base
SWI-Prolog and $BMLP_{active}$. Experimental phenotypic data on the same 16
carbon sources and FBA phenotypic predictions of iML1515 are already processed
as binary classes in iML1515 . In iML1515 , measured and predicted mutant
growth rates at a proportion of the wild-type’s growth rate greater than a
fixed threshold are mapped to the no phenotypic effect class.

Methods. We randomly sampled 100 simulations to perform and compute the
average wall time for each computation method. We first compared the mean and
standard deviation in simulation time of our new logical matrix encoding
against the Prolog encoding. In addition, we explored multi-threading with
both encodings. We used the multi-threading library offered by SWI-Prolog.
Owing to the overhead of managing threads, multi-threading boolean operations
with SWI-Prolog have diminishing returns. In BMLP-IE, simulation parameters
are encoded as matrices, enabling us to run simulations in batches. For BMLP-
IE, we used the concurrent threading library in SWI-Prolog at the batch level.
The runtime experiments with the multi-threading are repeated 10 times. We
show this comparison in Table 1.

Figure 4: Normalised confusion matrices of (a) BMLP-IE and (b) FBA phenotypic
predictions based on the GEM model iML1515. Each number has been normalised by
the size of the experiment dataset in iML1515 . From this dataset, the
phenotypic prediction accuracies of BMLP-IE and FBA are $88.7\%$ and $93.4\%$
respectively.

To evaluate BMLP-IE’s accuracy, we compared its phenotypic predictions against
experimental phenotypic data on the same 16 carbon sources. We identified the
number of consistent and inconsistent predictions and created normalised
confusion matrices for BMLP-IE and the FBA predictions from iML1515 . We show
these normalised confusion matrices in Fig 4.

Results and discussion. Table 1 shows a 170 times improvement in simulation
time by BMLP-IE compared to base SWI-Prolog. The simulation runtime is further
improved with multi-threading. With 20 cores, BMLP-IE uses only
$\frac{1}{600}$ of the simulation time compared to running simulations with
base SWI-Prolog. However, the multi-threading runtime improvement is not as
pronounced as the runtime improvement from using BMLP-IE. This result answers
Q1 as the SWI-Prolog implementation of BMLP-IE significantly improves the
phenotype simulation runtime of iML1515. While BMLP-IE is a specifical-purpose
algorithm, this result shows that a machine learning system implemented in
SWI-Prolog can be more time-efficient if it uses BMLP-IE when the background
knowledge can be described as a linear and immediately recursive datalog
program.

In Fig 4, the confusion matrices show a comparable predictive accuracy by
BMLP-IE ($88.7\%$) to that of standard FBA ($93.4\%$). This answers the
question Q2. Compared with FBA’s growth rate predictions, the
misclassifications by BMLP-IE are mainly due to it not being able to predict
reduced cell growth. FBA permit quantitative descriptions of metabolic fluxes
but metabolites are only binarily represented in an OEN. However, the strength
of the BMLP-IE approach is that it is easily adaptable for extending the
underlying GEM. FBA methods are mainly used as mathematical modelling tools so
they are not immediately applicable for learning.

###  8.2 Experiment 2: active learning sample complexity and experiment cost

Materials. Up to 3 of 16 common carbon sources were used as additional
nutrients to the base growth medium. Each of these 16 carbon sources is an
optional nutrient in the experimental data for validating iML1515 . This set
contains preferred carbon sources martinez-gomez_new_2012 ;
zampieri_regulatory_2019 for E. coli such as glucose and glycerol. We first
simulated iML1515 with up to 3 optional carbon sources, generating
$\sum_{0}^{i=3}\binom{16}{i}=697$ synthetic data. This allows us to go beyond
the 16 experimental data for each gene knockout in iML1515 . We observed a
single-knockout phenotypic effect from 213 genes. We treated these synthetic
data as ground truths where 213 genes have only positive examples. We only
focused on a subset of 2 genes because they have both positive and negative
examples555The distribution of training examples for this experiment is
described in the Appendix.. Then, we removed gene functions associated with
these 2 genes from the background knowledge. Thus, the resulting incomplete
model has a smaller set of $codes/2$ clauses.

Methods. We focused on one removed gene function at a time to recover it. Each
hypothesis except the empty hypothesis associates a gene locus identifier,
whose function we tried to recover, with an $enzyme/6$ clause. We concentrate
on 1294 $enzyme/6$ clauses that represent the non-maintenance reactions in
iML1515. Therefore, the hypothesis space had 1295 $codes/2$ hypotheses,
including an empty hypothesis. The instance space contained 697 experiments
with optional growth nutrients up to 3 carbon sources.

Figure 5: Total experimentation cost needed for the recovered models to reach
different levels of predictive accuracy. Models recovered from random
experiment selection and $BMLP_{active}$ can both converge to 100% predictive
accuracy. However, models recovered by $BMLP_{active}$ reach 100% predictive
accuracy with $\frac{1}{10}$ reagent cost expenditure of random experiment
sampling. The predictive accuracy of an empty hypothesis is 17.4%.

$BMLP_{active}$ actively selected experiments from the instance space to
reduce the hypothesis space. The random selection method randomly sampled
instances from the instance space. First, we constrained the selection of
experiments by an experimental reagent cost budget $C$. After experiment
selection, both methods returned the hypothesis with the single highest
compression, otherwise, a hypothesis was selected from those sharing the
highest compression randomly. The final hypotheses obtained by either method
would output a recovered model. We compared $BMLP_{active}$ with the random
experiment selection strategy with increasing experimental reagent cost budget
$C\in\\{1,10,50,100,500,1000,5000,10000\\}$. We tested the predictive accuracy
of the recovered models against all 697 synthetic data. The cost of selected
experiments and predictive accuracy of the recovered model were recorded for
$BMLP_{active}$ and random selection. Each selection method was repeated 10
times for each reagent cost budget $C$. We show this comparison in Fig 5.

Figure 6: Reduction of the number of experiments needed to recover gene
functions (sample complexity reduction). $BMLP_{active}$ learns 100% accurate
gene functions with 3 experiments while random sampling requires 25
experiments to do so. The predictive accuracy of an empty hypothesis is 17.4
%. Predictive accuracy is higher for both methods initially since experiment
selection is no longer constrained by reagent costs.

In addition, we constrained the number of experiments $N$ that can be selected
by increasing it linearly according to $\\{0,1,2,3,4,5,10,15,20,25,30\\}$. We
applied $BMLP_{active}$ and random experiment selection with compression to
decide the final hypothesis for recovering models. Each selection method was
repeated 10 times for every number of selected experiments $N$. Then, we
computed the phenotypic predictive accuracy of the recovered models against
all 697 synthetic data. We show this comparison in Fig 6. For experiment
number limits excluding 0, we calculated the mean error of learned hypotheses
and divided it by the number of experiments selected to obtain the ratio. In
Fig 7, we demonstrate how these ratios change to support Theorem 6.1.

Results and discussion. Fig 5 shows that $BMLP_{active}$ spends 10% of the
cost used by random experiment sampling when converging to the correct
hypothesis. The y-axis is the total experimental reagent cost in $log_{10}$
calculated from selected experiments. And 100% predictive accuracy indicates
successful recovery of the deleted gene functions. The correct gene functions
were recovered by $BMLP_{active}$ and random experiment selection with
$10^{1.99}$ and $10^{2.99}$ experimental reagent costs respectively. The total
experimental reagent cost is computed by dividing the unit cost which is the
cost of the cheapest optional nutrient ($\pounds 0.038/gram$ ). Given the
current hypothesis space and experimental design space and assuming each
reagent was used one gram, experiments selected by $BMLP_{active}$ and random
sampling would cost $10^{1.99}\times 0.038\approx\pounds 3.8$ and
$10^{2.99}\times 0.038\approx\pounds 38$ to learn a single gene function.
Consideration of differential resource factors in the user-defined
experimental cost function could significantly increase the total experimental
cost.

We then examined the advantage of $BMLP_{active}$ over random experiment
selection on the number of experiments needed to recover gene functions. The
number of experiments is a significant determinant of experimental cost. Fig 6
shows that $BMLP_{active}$ reduces the number of instances needed for
labelling to learn accurate gene functions. While random experiment selection
requires 25 experiments to recover the deleted gene functions, $BMLP_{active}$
only needed 3 experiments. Fig 5 and Fig 6 answer the question Q3 since
$BMLP_{active}$ successfully reduced the experimental reagent cost and
experiment cost compared to random experiment selection. Fig 5 and Fig 6 show
$BMLP_{active}$ can simultaneously optimise a user-defined experimental cost
function and reduce the total number of experiments to learn gene functions
accurately. This means that $BMLP_{active}$ is highly flexible and could be
tailored for specific experimental objectives in a discovery process. Later,
in Experiment 3, we revisit the question Q3 for recovering a known gene-
isoenzyme association which presents a bigger empirical challenge due to the
need to explore a larger experimental design space.

Figure 7: Ratio of sample complexity. $BMLP_{active}$ selects instances with a
minimal reduction ratio of at least 0.5%. For error $\epsilon>0$, the ratio of
$BMLP_{active}$ is clearly above the ratio of random sampling. This shows that
Theorem 6.1 correctly describes the sample complexity relationship between
active learning and random sampling.

Fig 7 demonstrates the mathematical relationships between predictive error and
the number of experiments required by $BMLP_{active}$ and random experiment
selection to converge to the correct hypothesis. For this gene function
recovery task, $BMLP_{active}$ selected experiments with minimal reduction
ratio $\phi\geq 0.5\%$. The ratio $\frac{\epsilon}{s_{active}}$ is greater
than $\frac{\epsilon+\phi}{s_{passive}}$ for error $\epsilon>0$. There are
fewer points on the active learning ratio curve since $BMLP_{active}$
converges to 100% accuracy using fewer training examples. The two curves
intersect at $\epsilon=0$ since both converged to 100% accuracy. Fig 7 answers
the question Q4 and supports Theorem 6.1. Theorem 6.1 shows that active
learning ($BMLP_{active}$) can reduce the requirement for training examples
compared with passive learning (random experiment selection) to achieve a
certain learning performance. In the context of scientific discovery, the
experimental cost to arrive at a finding can be reduced by active learning
since each actively selected experiment can have a higher return of
information.

###  8.3 Experiment 3: learning digenic functions

Materials. We focused on re-discovering a target isoenzyme function associated
with the gene $tyrB$ in the Tryptophan biosynthesis pathways. We aim to
demonstrate gene-isoenzyme function learning by recovering this gene function.
Associated functions of $tyrB$ and $aspC$ can be classified as isoenzyme
functions since they are responsible for producing an important amino acid
L-phenylalanine. In our logical representation of iML1515, we use reaction and
$codes/2$ unit ground literals to describe an isoenzyme. We selected 33 genes
associated with this essential pathway. In this focus gene set, 27 genes are
associated with a single function and 6 genes are related to two functions. We
consider medium conditions with 7 optional nutrients including aromatic amino
acids for recovering the deleted $tyrB$ isoenzyme function. Aromatic amino
acids are fundamental for protein synthesis but also serve as precursors for
vital secondary metabolites and high-value chemicals choi_violacein_2015 .

Methods. We removed tyrB’s gene function and metabolic reaction associated
with this isoenzyme from $BMLP_{active}$. The experiment instance space
contained double gene-knockout and single gene-knockout experiments of 33
genes. There were $\binom{33}{2}$ gene pair combinations for double-gene
knockout experiments and 33 single-gene knockout experiments. We considered
growth medium conditions with the 7 optional nutrients. Therefore, the
experiment instance space was $(\binom{33}{2}+33)\times 7=3927$ and 3696
experiment instance labels were synthetically generated from double-knockout
experiments using the full iML1515 background knowledge. The rest 231 single-
knockout empirical data were from validation data in iML1515 .

In the metabolic network, an isoenzyme function hypothesis contains a node
represented by a new enzyme identifier (a $enzyme/6$ clause) and a new link
between a gene locus identifier and this new enzyme identifier (a $codes/2$
clause). The hypothesis space contained associations of 33 genes (from the
focused library) with 33 $enzymes/6$ clauses. There were $(27\times 32+6\times
31)=1050$ potential hypotheses for 27 single-function genes and 6 double-
function genes. We also included the tyrB function hypothesis and an empty
hypothesis, which gave 1052 hypotheses in total. We increased linearly and
exponentially the number of experiments $N$ that could be selected according
to {1, 2, 4, 5, 8, 10, 15, 16, 20, 25, 30, 32, 64, 128, 256}. $BMLP_{active}$
selected $N$ experiments from this instance space to actively learn from the
hypothesis space. The random selection strategy randomly sampled $N$ instances
from the instance space. Again, both methods output the hypothesis with the
highest compression according to the sampled examples. A hypothesis was
randomly selected when multiple competing hypotheses shared the highest
compression. We recovered the models by adding the final hypotheses. We
recorded the number of experiments selected and evaluated if the recovered
models contained the target isoenzyme function from the iML1515 background
knowledge. Each experiment selection method was repeated 10 times and we
computed the frequency of successful isoenzyme recovery across these 10
repeats. We show this comparison in Fig 8.

Figure 8: tyrB isoenzyme function recovery frequency. We obtain the frequency
of successful recovery from various randomisation seeds and compute the
standard errors. We observe that within 64 experiments, random experiment
selection does not recover the tyrB isoenzyme function. In contrast,
$BMLP_{active}$ can recover this function with high frequency with 20
experiments.

Results and discussion. Fig 8 shows the result of recovering tyrB’s isoenzyme
function. For answering question Q5, we observed that $BMLP_{active}$
successfully recovers the correct tyrB isoenzyme function in 20 experiments.
This demonstrates $BMLP_{active}$’s ability to discover gene-isoenzyme
associations in GEM. However, to learn all isoenzyme functions, a high-
throughput experimentation procedure is required and we further discuss this
in Section 9.

This result also shows that $BMLP_{active}$ significantly reduced the number
of experiments needed compared to random experiment selection, which answers
the question Q3. $BMLP_{active}$ provides higher information gain from each
experiment and can guarantee recovery with as few as 20 experiments. While
random sampling could prune most candidate hypotheses, it failed to eliminate
competitive hypotheses further. In contrast to $BMLP_{active}$, random
sampling could not recover this isoenzyme function with 64 experiments and
could only do this occasionally with more than 10 times the number of
experiments required by $BMLP_{active}$. Although we used a reduced set of 30
genes in this demonstration, the combinatorial space was sufficiently large
that random experiment selection became non-viable as an experimentation
strategy in a discovery process. The remarkable increase in efficiency with
$BMLP_{active}$ with this task demonstrates it has potential even as the
experimental design space grows.

##  9 Conclusion and future work

We present a new logic programming system $BMLP_{active}$ that can simulate
and learn gene functions in a state-of-the-art GEM. This GEM, iML1515,
describes 1515 genes and 2719 metabolic reactions in E. coli. To our
knowledge, $BMLP_{active}$ is the first logic programming system able to
operate on GEMs. We implement $BMLP_{active}$ via our proposed BMLP approach
by exploiting boolean matrices to improve the simulation runtime of a GEM by
170 times. Our empirical results show that $BMLP_{active}$ can learn gene
functions while reducing the cost of optional nutrient substance in
experiments by 90%. In addition, we deploy $BMLP_{active}$ to actively learn
gene-isoenzyme functions which are an error source in GEM. As a demonstration,
we empirically show that $BMLP_{active}$ can learn an important isoenzyme
function in the aromatic amino acid biosynthesis pathway with as little as 20
training examples. In comparison, random experiment selection cannot learn
this isoenzyme function consistently even with more than 250 training data.
$BMLP_{active}$ provides a substantial advantage for isoenzyme function
discovery over random experiment selection due to the sample complexity
reduction. This finding suggests that $BMLP_{active}$ could be used to offset
the growth in experimental design spaces associated with genome-scale function
learning.

One limitation of $BMLP_{active}$ is that it cannot yet quantitatively predict
the behaviours of biological systems. Quantitative growth predictions would
allow us to model in-vivo phenotype observations without binary
discretisation. Future work will focus on using logic programs to optimise
relations and parameters given numeric training data. A natural link exists
between Petri nets and Probabilistic Logic Programming
de_raedt_probabilistic_2015 . Transitions sometimes have firing likelihood
constraints and the uncertainty can be modelled by probabilistic logic
programs. Evaluating these programs requires tensor-based algebra of logic
programs, a topic that has been explored recently in sakama_logic_2021 .

$BMLP_{active}$ also requires a pre-defined hypothesis space since it does not
invent predicates for learning new concepts. This relates to a central problem
of ILP known as predicate invention ILP30 . New predicate names and reusable
definitions automatically learned in addition to the background knowledge are
machine inventions. Recent methods have combined ILP with meta-level logic
programs and matrices such as the Deeplog system muggleton_hypothesizing_2023
, $\delta$ILP evans_learning_2018 and $Meta_{Abd}$ dai_abductive_2021 .
Reusing invented predicates can reduce the size of a hypothesis learned by an
ILP system and improve its learning performance ILP30 . $BMLP_{active}$ could
be extended with predicate invention for learning novel biological knowledge.

To make a biological discovery, we aim to integrate this approach with an
experiment workflow such as the one shown in wang_pooled_2018 to address
whole-genome function learning. Data created from a high-throughput empirical
framework can be first compared with simulated phenotypic predictions. The
inconsistency can inform us of the error sources in the biological model,
narrowing the scope of hypothesis space. With advanced biological gene editing
tools, genes edition can be targeted for experimental confirmation of
candidate hypotheses during active learning. Importantly, we show that our
approach can address the interactions between multiple genes, which are not
addressable by current mutagenesis approaches.

Designing experiments to explain some hypothesised phenomenon requires
deliberate thinking by scientists. A helpful scientific assistant could teach
scientists to create experiment plans with optimal outcomes. This aspect of
human-AI collaboration has been studied in an emerging research area on Ultra-
Strong Machine Learning (USML) michie1988 . A USML not only learns a concept
but also teaches this concept to humans such that human comprehension is
improved. Prior research has shown instances of USML in the context of
learning recursive concepts US2018 and a two-player game Ai2021 . A recent
study sequential_teaching demonstrated that an ILP system can play the role of
a “teacher” and human learning from machine explanations of its output can
sometimes inspire human discovery. To date, there is no demonstration of
autonomous discovery systems that satisfy the USML criterion. The ultimate
goal of autonomous discovery systems is to acquire new scientific
understanding and explain insights to humans krenn_scientific_2022 . The
present work can be extended to examine the effect of machine-selected
experiments and machine-learned output on human experiment selection and human
comprehension of complex biological concepts.

## Ethical Statement

Datasets used do not contain any personal information or other sensitive
content. There are no ethical issues.

## Code availability

Data used and generated in all experimental sections are accessible on
https://github.com/lai1997/BMLP_active_public.

## Acknowledgments

The first, third and fourth authors acknowledge support from the UKRI 21EBTA:
EB-AI Consortium for Bioengineered Cells & Systems (AI-4-EB) award
(BB/W013770/1). The second author acknowledges support from the UK’s EPSRC
Human-Like Computing Network (EP/R022291/1), for which he acts as Principal
Investigator. ChatGPT and Grammarly have been used as language editing tools
after all intellectual content has been drafted.

## References

  * [1] L. Ai, J. Langer, S. H. Muggleton, and U. Schmid.  Explanatory machine learning for sequential human teaching.  Machine Learning, 112:3591–3632, 2023.
  * [2] L. Ai, S. Muggleton, C. Hocquette, M. Gromowski, and U. Schmid.  Beneficial and harmful explanatory machine learning.  Machine Learning, 110:695–721, 2021.
  * [3] C. Angione.  Human Systems Biology and Metabolic Modelling: A Review—From Disease Metabolism to Precision Medicine.  BioMed Research International, 2019:e8304260, 2019.
  * [4] D. Angluin.  Queries Revisited.  In Proceedings of the 12th International Conference on Algorithmic Learning Theory, ALT ’01, pages 12–31, Berlin, Heidelberg, 2001. Springer-Verlag.
  * [5] S. Anwar, C. Baral, and K. Inoue.  Encoding Higher Level Extensions of Petri Nets in Answer Set Programming.  In P. Cabalar and T. C. Son, editors, Logic Programming and Nonmonotonic Reasoning, Lecture Notes in Computer Science, pages 116–121, Berlin, Heidelberg, 2013. Springer.
  * [6] G. W. Beadle and E. L. Tatum.  Genetic Control of Biochemical Reactions in Neurospora.  Proceedings of the National Academy of Sciences, 27(11):499–506, 1941.
  * [7] T. M. Behrens and J. Dix.  Model checking multi-agent systems with logic based Petri nets.  Annals of Mathematics and Artificial Intelligence, 51(2):81–121, 2007.
  * [8] C. Berge.  Hypergraphs: combinatorics of finite sets.  North-Holland mathematical library. North Holland Distributors for the U.S.A. and Canada, Elsevier Science Pub. Co, Amsterdam New York, 1989.
  * [9] D. B. Bernstein, B. Akkas, M. N. Price, and A. P. Arkin.  Evaluating E. coli genome-scale metabolic model accuracy with high-throughput mutant fitness data.  Molecular Systems Biology, 19(12), 2023.
  * [10] A. Blumer, A. Ehrenfeucht, D. Haussler, and M. K. Warmuth.  Occam’s Razor.  Information Processing Letters, 24(6):377–380, 1987.
  * [11] D. A. Boiko, R. MacKnight, B. Kline, and G. Gomes.  Autonomous chemical research with large language models.  Nature, 624(7992):570–578, 2023.
  * [12] A. M. Bran, S. Cox, O. Schilter, C. Baldassari, A. White, and P. Schwaller.  Augmenting large language models with chemistry tools.  In NeurIPS 2023 AI for Science Workshop, 2023.
  * [13] S. L. Brunton, J. L. Proctor, and J. N. Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, 2016.
  * [14] C. H. Bryant, S. H. Muggleton, S. G. Oliver, D. Kell, P. Reiser, and R. D. King.  Combining inductive logic programming, active learning and robotics to discover the function of genes.  Electronic Transactions in Artificial Intelligence, pages 1–36, 2001.
  * [15] S. Ceri, G. Gottlob, and L. Tanca.  What you always wanted to know about Datalog (and never dared to ask).  IEEE Transactions on Knowledge and Data Engineering, 1(1):146–166, 1989.
  * [16] S. Y. Choi, K.-h. Yoon, J. I. Lee, and R. J. Mitchell.  Violacein: Properties and Production of a Versatile Bacterial Pigment.  BioMed Research International, 2015:1–8, 2015.
  * [17] W. Cohen, F. Yang, and K. R. Mazaitis.  TensorLog: A Probabilistic Database Implemented Using Deep-Learning Infrastructure.  Journal of Artificial Intelligence Research, 67:285–325, 2020.
  * [18] D. Cohn, L. Atlas, and R. Ladner.  Improving generalization with active learning.  Machine Learning, 15(2):201–221, 1994.
  * [19] D. Conklin and I. H. Witten.  Complexity-based induction.  Machine Learning, 16(3):203–225, 1994.
  * [20] I. M. Copilowish.  Matrix Development of the Calculus of Relations.  The Journal of Symbolic Logic, 13(4):193–203, 1948.
  * [21] M. Costanzo, E. Kuzmin, J. van Leeuwen, B. Mair, J. Moffat, C. Boone, and B. Andrews.  Global Genetic Networks and the Genotype-to-Phenotype Relationship.  Cell, 177(1):85–100, 2019.
  * [22] M. Costanzo, B. VanderSluis, E. N. Koch, et al.  A global genetic interaction network maps a wiring diagram of cellular function.  Science, 353(6306), 2016.
  * [23] A. Cropper, S. Dumancic, R. Evans, and S. H. Muggleton.  Inductive logic programming at 30.  Machine Learning, 111:147–172, 2021.
  * [24] W.-Z. Dai and S. H. Muggleton.  Abductive Knowledge Induction from Raw Data.  In Z.-H. Zhou, editor, Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21, pages 1845–1851. International Joint Conferences on Artificial Intelligence Organization, 2021.
  * [25] S. Dasgupta.  Coarse sample complexity bounds for active learning.  In Advances in Neural Information Processing Systems, volume 18. MIT Press, 2005.
  * [26] L. De Raedt and A. Kimmig.  Probabilistic (logic) programming concepts.  Machine Learning, 100(1):5–47, 2015.
  * [27] Y. Dimopoulos, E. Kouppari, A. Philippou, and K. Psara.  Encoding Reversing Petri Nets in Answer Set Programming.  In I. Lanese and M. Rawski, editors, Reversible Computation, Lecture Notes in Computer Science, pages 264–271, Cham, 2020. Springer International Publishing.
  * [28] A. Domenici.  Petri nets in logic.  Microprocessing and Microprogramming, 30(1):193–198, 1990.
  * [29] R. Evans and E. Grefenstette.  Learning Explanatory Rules from Noisy Data.  Journal of Artificial Intelligence Research, 61:1–64, 2018.
  * [30] L. Faure, B. Mollet, W. Liebermeister, and J.-L. Faulon.  A neural-mechanistic hybrid approach improving the predictive power of genome-scale metabolic models.  Nature Communications, 14(1):4669, 2023.
  * [31] M. J. Fischer and A. R. Meyer.  Boolean matrix multiplication and transitive closure.  In 12th Annual Symposium on Switching and Automata Theory (SWAT 1971), pages 129–131, 1971.
  * [32] M. Gelfond and V. Lifschitz.  The Stable Model Semantics for Logic Programming.  In Proceedings of International Logic Programming Conference and Symposium, pages 1070–1080, 1988.
  * [33] E. Grefenstette.  Towards a Formal Distributional Semantics: Simulating Logical Calculi with Tensors.  In Proceedings of the Second Joint Conference on Lexical and Computational Semantics (SEM), pages 1–10, 2013.
  * [34] R. Guimerà, I. Reichardt, A. Aguilar-Mogas, F. A. Massucci, M. Miranda, J. Pallarès, and M. Sales-Pardo.  A Bayesian machine scientist to aid in the solution of challenging scientific problems.  Science Advances, 6(5), 2020.
  * [35] D. Heckmann, C. J. Lloyd, N. Mih, Y. Ha, D. C. Zielinski, Z. B. Haiman, A. A. Desouki, M. J. Lercher, and B. O. Palsson.  Machine learning applied to enzyme turnover numbers reveals protein structural correlates and improves metabolic models.  Nature Communications, 9(1):5252, 2018.
  * [36] C. Hocquette and S. Muggleton.  How Much Can Experimental Cost Be Reduced in Active Learning of Agent Strategies?  In F. Riguzzi, E. Bellodi, and R. Zese, editors, Inductive Logic Programming, volume 11105, pages 38–53. Springer International Publishing, Cham, 2018.
  * [37] Y. E. Ioannidis.  On the Computation of the Transitive Closure of Relational Operators.  In Proceedings of the 12th International Conference on Very Large Data Bases, VLDB ’86, pages 403–411, San Francisco, CA, USA, 1986.
  * [38] R. D. King, J. Rowland, S. G. Oliver, M. Young, W. Aubrey, E. Byrne, M. Liakata, M. Markham, P. Pir, L. N. Soldatova, A. Sparkes, K. E. Whelan, and A. Clare.  The Automation of Science.  Science, 324(5923):85–89, 2009.
  * [39] R. D. King, K. E. Whelan, F. M. Jones, P. G. K. Reiser, C. H. Bryant, S. H. Muggleton, D. B. Kell, and S. G. Oliver.  Functional genomic hypothesis generation and experimentation by a robot scientist.  Nature, 427:247–252, 2004.
  * [40] Z. A. King, J. Lu, A. Dräger, P. Miller, S. Federowicz, J. A. Lerman, A. Ebrahim, B. O. Palsson, and N. E. Lewis.  BiGG Models: A platform for integrating, standardizing and sharing genome-scale models.  Nucleic Acids Research, 44(D1):D515–D522, 2016.
  * [41] M. Krenn, R. Pollice, S. Y. Guo, and others.  On scientific understanding with artificial intelligence.  Nature Reviews Physics, 4:761–769, 2022.
  * [42] P. W. Langley, H. A. Simon, G. Bradshaw, and J. M. Zytkow.  Scientific Discovery: Computational Explorations of the Creative Process.  The MIT Press, 1987.
  * [43] F. Lin.  From Satisfiability to Linear Algebra.  Technical report, Hong Kong University of Science and Technology, 2013\.
  * [44] J. W. Lloyd.  Foundations of Logic Programming.  Springer Science and Business Media, 2012.
  * [45] K. Martínez-Gómez, N. Flores, H. M. Castañeda, G. Martínez-Batallar, G. Hernández-Chávez, O. T. Ramírez, G. Gosset, S. Encarnación, and F. Bolivar.  New insights into Escherichia coli metabolism: carbon scavenging, acetate metabolism and carbon recycling responses during growth on glycerol.  Microbial Cell Factories, 11(1):46, 2012.
  * [46] D. Michie.  Machine learning in the next five years.  In Proceedings of the 3rd European Working Session on Learning, pages 107–122. Pitman, 1988.
  * [47] T. M. Mitchell.  Generalization as search.  Artificial Intelligence, 18:203–226, 1982.
  * [48] J. M. Monk, C. J. Lloyd, E. Brunk, N. Mih, A. Sastry, Z. King, R. Takeuchi, W. Nomura, Z. Zhang, H. Mori, A. M. Feist, and B. O. Palsson.  iML1515, a knowledgebase that computes escherichia coli traits.  Nature Biotechnology, 35(10):904–908, 2017.
  * [49] S. H. Muggleton.  Inductive logic programming.  New Generation Computing, 8:295–318, 1991.
  * [50] S. H. Muggleton.  Inverse entailment and Progol.  New Generation Computing, 13:245–286, 1995.
  * [51] S. H. Muggleton.  Hypothesizing an algorithm from one example: the role of specificity.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 381(2251):20220046, 2023.
  * [52] S. H. Muggleton and C. H. Bryant.  Theory completion using inverse entailment.  In Proceedings of the 10th International Workshop on Inductive Logic Programming (ILP-00), pages 130–146, 2000.
  * [53] S. H. Muggleton, U. Schmid, C. Zeller, A. Tamaddoni-Nezhad, and T. Besold.  Ultra-strong machine learning: Comprehensibility of programs learned with ILP.  Machine Learning, 107:1119–1140, 2018.
  * [54] S. Nienhuys-Cheng and R. Wolf.  Foundations of Inductive Logic Programming. Springer-Verlag New York, Inc., 1997.
  * [55] T. Oyetunde, D. Liu, H. G. Martin, and Y. J. Tang.  Machine learning framework for assessment of microbial factory performance.  PLOS ONE, 14(1):e0210558, 2019.
  * [56] B. O. Palsson.  Systems Biology: Constraint-based Reconstruction and Analysis.  Cambridge University Press, Cambridge, 2015.
  * [57] C. S. Peirce.  Collected Papers of Charles Sanders Peirce, Volumes II.  Harvard University Press, 1932.
  * [58] B. K. Petersen, M. L. Larma, T. N. Mundhenk, C. P. Santiago, S. K. Kim, and J. T. Kim.  Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients.  In Proceedings of the 9th International Conference on Learning Representations, 2020.
  * [59] M. N. Price, G. M. Zane, J. V. Kuehl, R. A. Melnyk, J. D. Wall, A. M. Deutschbauer, and A. P. Arkin.  Filling gaps in bacterial amino acid biosynthesis pathways with high-throughput genetics.  PLOS Genetics, 14(1):e1007147, 2018.
  * [60] P. Rana, C. Berry, P. Ghosh, and S. S. Fong.  Recent advances on constraint-based models by integrating machine learning.  Current Opinion in Biotechnology, 64:85–91, 2020.
  * [61] W. Reisig.  Understanding Petri Nets: Modeling Techniques, Analysis Methods, Case Studies.  Springer, Berlin, Heidelberg, 2013.
  * [62] Z. Ren, Z. Ren, Z. Zhang, T. Buonassisi, and J. Li.  Autonomous experiments using active learning and AI.  Nature Reviews Materials, 8(9):563–564, 2023.
  * [63] G. Rozenberg and J. Engelfriet.  Elementary net systems.  In Reisig, W., Rozenberg, G. (eds) Lectures on Petri Nets I: Basic Models. ACPN 1996. Lecture Notes in Computer Science, volume 1491. Springer Berlin Heidelberg, 1998.
  * [64] A. Sahu, M.-A. Blätke, J. J. Szymański, and N. Töpfer.  Advances in flux balance analysis by integrating machine learning and mechanism-based models.  Computational and Structural Biotechnology Journal, 19:4626–4640, 2021.
  * [65] C. Sakama, K. Inoue, and T. Sato.  Logic programming in tensor spaces.  Annals of Mathematics and Artificial Intelligence, 89(12):1133–1153, 2021.
  * [66] T. Sato.  A linear algebraic approach to datalog evaluation.  Theory and Practice of Logic Programming, 17(3):244–265, 2017.
  * [67] T. Sato and K. Inoue.  Differentiable learning of matricized DNFs and its application to Boolean networks.  Machine Learning, 112(8):2821–2843, 2023.
  * [68] T. Sato, K. Inoue, and C. Sakama.  Abducing Relations in Continuous Spaces.  In Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, pages 1956–1962, 2018.
  * [69] T. Sato and R. Kojima.  Boolean Network Learning in Vector Spaces for Genome-wide Network Analysis.  In Proceedings of the Eighteenth International Conference on Principles of Knowledge Representation and Reasoning, pages 560–569, Hanoii, Vietnam, 2021. International Joint Conferences on Artificial Intelligence Organization.
  * [70] P. Sen, S. Lamichhane, V. B. Mathema, A. McGlinchey, A. M. Dickens, S. Khoomrung, and M. Orešič.  Deep learning meets metabolomics: a methodological perspective.  Briefings in Bioinformatics, 22(2):1531–1542, 2021.
  * [71] C. Shannon and W. Weaver.  The Mathematical Theory of Communication.  University of Illinois Press, Urbana, 1963.
  * [72] A. Srinivasan.  The ALEPH manual.  Machine Learning at the Computing Laboratory, Oxford University, 2001.
  * [73] A. Srinivasan and M. Bain.  Knowledge-Guided Identification of Petri Net Models of Large Biological Systems.  In S. H. Muggleton, A. Tamaddoni-Nezhad, and F. A. Lisi, editors, Inductive Logic Programming, Lecture Notes in Computer Science, pages 317–331, 2012.
  * [74] A. Tarski.  A lattice-theoretical fixpoint theorem and its applications.  Pacific Journal of Mathematics, 5(2):85–309, 1955.
  * [75] L. Todorovski and S. Džeroski.  Integrating knowledge-driven and data-driven approaches to modeling.  Ecological Modelling, 194(1):3–13, 2006.
  * [76] C. Tosh and S. Dasgupta.  Diameter-Based Active Learning.  In Proceedings of the 34th International Conference on Machine Learning, pages 3444–3452. PMLR, 2017.
  * [77] M. H. Van Emden and R. A. Kowalski.  The Semantics of Predicate Logic as a Programming Language.  Journal of the ACM, 23(4):733–742, 1976.
  * [78] T. Wang, C. Guan, J. Guo, B. Liu, Y. Wu, Z. Xie, C. Zhang, and X.-H. Xing.  Pooled CRISPR interference screening enables genome-scale functional genomics study in bacteria with superior performance.  Nature Communications, 9(1):2475, 2018.
  * [79] S. G. Wu, Y. Wang, W. Jiang, T. Oyetunde, R. Yao, X. Zhang, K. Shimizu, Y. J. Tang, and F. S. Bao.  Rapid Prediction of Bacterial Heterotrophic Fluxomics Using Machine Learning and Constraint Programming.  PLOS Computational Biology, 12(4), 2016.
  * [80] F. Yang, Z. Yang, and W. W. Cohen.  Differentiable learning of logical rules for knowledge base reasoning.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, pages 2316–2325, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [81] J. H. Yang, S. N. Wright, M. Hamblin, D. McCloskey, M. A. Alcantar, L. Schrübbers, A. J. Lopatkin, S. Satish, A. Nili, B. O. Palsson, G. C. Walker, and J. J. Collins.  A White-Box Machine Learning Approach for Revealing Antibiotic Mechanisms of Action.  Cell, 177(6):1649–1661, 2019.
  * [82] M. Zampieri, M. Hörl, F. Hotz, N. F. Müller, and U. Sauer.  Regulatory mechanisms underlying coordination of amino acid and glucose catabolism in Escherichia coli.  Nature Communications, 10(1):3354, 2019.

## Appendix A. Reagent (optional nutrient substance) costs

Carbon Source | Cost per gram (£/g)
---|---
D-Mannitol | 0.534
Xylose | 0.51
Galactose | 0.0936
Pyruvate | 0.2776
Ribose | 0.993
Oxoglutarate | 0.437
Acetate | 0.0748
Gluconate | 0.0359
Glycerol | 0.07
Maltose | 0.244
Lactate | 0.0546
Sorbitol | 0.196
Succinate | 0.0825
N-acetyl-glucosamine | 2.93
L-Alanine | 0.81
Glucose | 0.0476
Table 2: The costs of 16 optional carbon sources considered in Experiment 2. Carbon Source | Cost per gram (£/g)
---|---
Indole | 0.0428
Phenylalanine | 0.0346
Tryptophan | 0.0346
Tyrosine | 0.0346
Shikimate | 8.2
L-Aspartic acid | 0.0156
Table 3: The costs of aromatic amino acids considered in Experiment 3.

## Appendix B. Gene function learning example distribution in Experiment 2.

Table 4: The distribution of synthetic positive and negative examples for each gene knockout based on simulation of iML1515. Gene identifiers start with the lower-case letter “b” and gene descriptive names are given in parentheses.  Gene Id. | $|E^{+}|$ | $|E^{-}|$
---|---|---
b0720 (gltA) | 576 | 121
b1136 (icd) | 576 | 121
b3729 (glmS) | 576 | 121
The other 210 genes | 697 | 0

In Experiment 2, for the 16 carbon sources, 213 genes have a single-knockout
phenotypic effect. Three genes, b0720 (gltA), b1136 (icd) and b3729 (glmS),
have both positive and negative examples. We focus on b0720 and b3729 since
b0720 and b1136 have the same coverage over the 697 training examples so their
functions cannot be distinguished.