# A Model Checking-based Analysis Framework for Systems Biology Models

Bing Liu  Department of Computational & Systems Biology
School of Medicine, University of Pittsburgh
Pittsburgh, PA 15237, U.S.A.
<EMAIL_ADDRESS>
###### Abstract

Biological systems are often modeled as a system of ordinary differential
equations (ODEs) with time-invariant parameters. However, cell signaling
events or pharmacological interventions may alter the cellular state and
induce multi-mode dynamics of the system. Such systems are naturally modeled
as hybrid automata, which possess multiple operational modes with specific
nonlinear dynamics in each mode. In this paper we introduce a model checking-
enabled framework than can model and analyze both single- and multi-mode
biological systems. We tackle the central problem in systems biology–identify
parameter values such that a model satisfies desired behaviors–using bounded
model checking. We resort to the delta-decision procedures to solve
satisfiability modulo theories (SMT) problems and sidestep undecidability of
reachability problems. Our framework enables several analysis tasks including
model calibration and falsification, therapeutic strategy identification, and
Lyapunov stability analysis. We demonstrate the applicablitliy of these
methods using case studies of prostate cancer progression, cardiac cell action
potential and radiation diseases.

###### Index Terms:

systems biology, model checking, hybrid systems, delta-decision, parameter
synthesis

##  I Introduction

Biomolecules interact with each other and form large and complicated networks.
A systems-wide view of functional regulation in the context of biochemical
networks is required for contemporary drug discovery and systems pharmacology
to identify new drug targets, understand drug adverse effects, and design
therapeutic strategies [1]. It is commonly recognized that systems modeling
will play a crucial role in this endeavor.

A standard approach of modeling the dynamics of a biochemical network is
through a system of ordinary differential equations (ODEs) [2]. However, model
parameters including kinetic rate constants and initial concentrations of
molecular species are often unknown and they will have to be estimated using
limited and noisy experimental data. Hence constructing and calibrating ODE
models of realistic biochemical networks is a challenging problem. Further,
ODE models that were developed for simulating single-cell behaviors were often
calibrated and validated using population-based data. This constitutes a
significant additional challenge.

To address these challenges, we have developed several model checking-enabled
analysis techniques for ODE models. For example, by reducing the ODE dynamics
as a dynamic Bayesian network (DBN) [3, 4, 5, 6], one can efficiently perform
parameter estimation and probabilistic model checking analysis by exploiting
Bayesian inferencing [7, 8]. This framework has been applied to study
complement system [9] and apoptosis pathway [10].

Figure 1: Distinct cell death pathways induced by irradiation exposure.

We also developed statistical model checking (SMC) techniques to calibrate and
analyze ODE systems with probabilistic initial states [11, 12, 13]. In this
setting, bounded linear temporal logic is used to encode quantitative
behavioral constraints and qualitative properties of biochemical networks. By
equipping existing parameter search algorithms with a SMC-based evaluation
method, one can arrive both novel and efficient parameter estimation and
analysis methods. This framework has been generalized to deal with stochastic
rule-based models [14, 15] and hybrid automata [16] and applied to studies of
various biological systems such as innate immune system [17] and cell
death/survival pathways [18, 19, 20, 21].

However, in many settings, it is not fruitful to view the functioning of a
biological system in terms of a single entity. Rather, the system will possess
multiple operational modes with a specific signaling network being active in
each mode. For example, signaling responses to ionizing irradiation exposure
within individual cells follow distinct cell death pathways (Fig. 1). The
interconnectivity between these pathways brings up new challenges on
determining the sequence and timing of medication administration against
radiation injuries, since interventions may alter the cellular state and
induce differential modes of dynamics [22, 23, 24].

In such situations, using a monolithic approach will result in a messy and
very large model with time-invariant parameters. Consequently, the modeling
and analysis approach based on the notion of multi-mode biological networks
will be very useful. Multi-mode dynamics and their modeling appear frequently
in physical and engineering settings [25]. Multiple variants of the formalism
called hybrid automata [26] are often used to model biological processes [27,
28, 29, 30, 31, 32]. However such models are difficult to analyze and to get
around this, most efforts end up imposing severe restrictions on the dynamic
laws associated with the modes [33, 34, 35].

In this paper, we present a novel model checking-based framework for analyzing
single- or multi-mode biological systems with nonlinear dynamics (Fig. 2).
Given a dynamical system, we describe the set of states of interest as a
first-order logic formula and perform bounded model checking to determine
reachability of these states. We adapt an interval constrains propagation
(ICP) based algorithm to explore the parameter spaces and identify the sets of
parameters using which the model satisfies desired behavior constraints. Note
that determining the truth value of first-order sentences over the reals with
nonlinear real functions is a well-known undefinable problem. We use
$\delta$-decision procedures [36] to ask for answers that may have one-sided
$\delta$-bounded errors. If the model satisfies the desired behavior (e.g. the
model matches training and testing data), we can carry out analysis tasks such
as stability analysis or identify novel therapeutic strategies using
$\delta$-decision procedures. Otherwise, we will conduct SMC-based analysis
[12] to generate new hypotheses and refine the model structure iteratively. We
demonstrate the applicability of our framework by “proof-of-concept” studies
on prostate cancer and cardiac disorders [37, 38, 39].

Figure 2: A model checking-enabled modeling framework.

Turning to related work, a survey of modeling and analysis of biological
systems using hybrid models can be found in [40]. Formal verification of
hybrid systems is a well-established domain [41]. Analyzing the properties of
biochemical networks using model checking techniques is being actively pursued
by a number of groups [42, 43, 44, 45]. Of particular interest in our context
are parameter synthesis methods which identify parameter values for which some
qualitative behavior [46, 47, 48].

The rest of the paper is organized as follows. We briefly introduce model
checking in the next section. We describe the $\delta$-decision procedures in
Section III. We then present the applications of $\delta$-decision-based
analysis methods in Section IV. In the final section, we summarize and discuss
future work.

##  II Model Checking

Amir Pnueli introduced temporal logics into the world of program verification
[49]. The algorithmic verification procedure called model checking was
formulated by Clarke and Emerson and independently by Joseph Sifakis [50].
Briefly, the model checking procedure operates as follows. Given a model
$\mathcal{M}$ with initial state $s$, a model checker decides if a property
written as a temporal logic formula $\phi$ is satisfied, denoted as
$\mathcal{M},s\models\phi$. This can be done by: (i) constructing a finite
(state) transition system corresponding to $\mathcal{M}$ in which each state
represents a possible configuration of the system and each transition
represents an evolution of the system from one configuration to another, and
(ii) verify whether $\phi$ is satisfied by exhaustively exploring the set of
system executions. The key feature of the model checking procedure is that it
is fully automated. Further, if it is not the case that
$\mathcal{M},s\models\phi$ then the procedure will usually return as a
“counter-example” an execution due to which the property is not being met by
the system. This can serve as a powerful debugging tool for systems that are
large and complex. An excellent starting point for exploring this whole field
is [51].

##  III $\delta$-Decision Procedures

###  III-A $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-Formulas and
$\delta$-Decisions Over The Reals

In [36], we developed a theory of decision problems over the reals with
computable functions. Most common continuous real functions are computable,
including solutions of Lipschitz-continuous ODEs. In fact, the notion of
computability of real functions directly corresponds to whether they can be
numerically simulated. We write $\mathcal{F}$ to denote an arbitrary
collection of symbols representing computable functions over $\mathbb{R}^{n}$
for various $n$. We consider the first-order formulas with a signature
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}=\langle 0,1,\mathcal{F},>\rangle$.
Note that constants are seen as 0-ary functions in $\mathcal{F}$.
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formulas are evaluated in the
standard way over the corresponding structure
$\mathbb{R}_{\mathcal{F}}=\langle\mathbb{R},\mathcal{F},>\rangle$. We use
atomic formulas of the form $t(x_{1},...,x_{n})>0$ or $t(x_{1},...,x_{n})\geq
0$, where $t(x_{1},...,x_{n})$ are built up from functions in $\mathcal{F}$.
To avoid extra preprocessing of formulas, we give an explicit definition of
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formulas as follows.

######  Definition 1 ($\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-Formulas)

Let $\mathcal{F}$ be a collection of Type 2 functions, which contains at least
$0$, unary negation -, addition $+$, and absolute value $|\cdot|$. We define:

$\displaystyle t$ $\displaystyle:=x\;|\;f(t(\vec{x})),\mbox{ where
}f\in\mathcal{F}\mbox{, possibly constant};$ $\displaystyle\varphi$
$\displaystyle:=t(\vec{x})>0\;|\;t(\vec{x})\geq
0\;|\;\varphi\wedge\varphi\;|\;\varphi\vee\varphi\;|\;\exists
x_{i}\varphi\;|\;\forall x_{i}\varphi.$

In this setting $\neg\varphi$ is regarded as an inductively defined operation
which replaces atomic formulas $t>0$ with $-t\geq 0$, atomic formulas $t\geq
0$ with $-t>0$, swaps $\wedge$ and $\vee$, and swaps $\forall$ and $\exists$.
Implication $\varphi_{1}\rightarrow\varphi_{2}$ is defined as
$\neg\varphi_{1}\vee\varphi_{2}$.

######  Definition 2 (Bounded Quantifiers)

We define

$\displaystyle\exists^{[u,v]}x.\varphi$ $\displaystyle=_{df}\exists x.(u\leq
x\land x\leq v\wedge\varphi),$ $\displaystyle\forall^{[u,v]}x.\varphi$
$\displaystyle=_{df}\forall x.((u\leq x\land x\leq v)\rightarrow\varphi),$

where $u$ and $v$ denote $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$ terms whose
variables only contain free variables in $\varphi$, excluding $x$. It is easy
to check that
$\exists^{[u,v]}x.\varphi\leftrightarrow\neg\forall^{[u,v]}x.\neg\varphi$.

We say a sentence is bounded if it only involves bounded quantifiers.

######  Definition 3 (Bounded
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-Sentences)

A bounded $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-sentence is

$Q_{1}^{[u_{1},v_{1}]}x_{1}\cdots
Q_{n}^{[u_{n},v_{n}]}x_{n}\;\psi(x_{1},...,x_{n}).$

$Q_{i}^{[u_{i},v_{i}]}$s are bounded quantifiers, and $\psi(x_{1},...,x_{n})$
is a quantifier-free $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formula.

We write $\psi(x_{1},...,x_{n})$ as
$\psi[t_{1}(\vec{x})>0,...,t_{k}(\vec{x})>0;t_{k+1}(\vec{x})\geq
0,...,t_{m}(\vec{x})\geq 0]$ to emphasize that $\psi(\vec{x})$ is a Boolean
combination of the atomic formulas shown.

######  Definition 4 ($\delta$-Variants)

Let $\delta\in\mathbb{Q}^{+}\cup\\{0\\}$, and $\varphi$ an
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formula of the form

$\varphi:\ Q_{1}^{I_{1}}x_{1}\cdots
Q_{n}^{I_{n}}x_{n}\;\psi[t_{i}(\vec{x},\vec{y})>0;t_{j}(\vec{x},\vec{y})\geq
0],$

where $i\in\\{1,...k\\}$ and $j\in\\{k+1,...,m\\}$. The $\delta$-weakening
$\varphi^{\delta}$ of $\varphi$ is defined as the result of replacing each
atom $t_{i}>0$ by $t_{i}>-\delta$ and $t_{j}\geq 0$ by $t_{j}\geq-\delta$.
That is,

$\varphi^{\delta}:\ Q_{1}^{I_{1}}x_{1}\cdots
Q_{n}^{I_{n}}x_{n}\;\psi[t_{i}(\vec{x},\vec{y})>-\delta;t_{j}(\vec{x},\vec{y})\geq-\delta].$

We then have the following main decidability result.

######  Theorem 1 ($\delta$-Decidability)

Let $\delta\in\mathbb{Q}^{+}$ be arbitrary. There is an algorithm which, given
any bounded $\varphi$, correctly returns one of the following two answers:

  * •

$\phi$ is false ($\mathsf{unsat}$);

  * •

$\phi^{\delta}$ is true ($\delta$-$\mathsf{sat}$).

Note when the two cases overlap, either answer is correct.

We call this new decision problem the $\delta$-decision problem for
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-sentences.

######  Definition 5 ($\delta$-Complete Decision Procedures)

If an algorithm solves the $\delta$-decision problem correctly for a set $S$
of $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-sentences, we say it is
$\delta$-complete for $S$.

From $\delta$-decidability, $\delta$-complete decision procedures always exist
for bounded $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formulas. In practice, we
have shown that the combination of the DPLL(T) framework and ICP indeed gives
us a $\delta$-complete decision procedure [52].

###  III-B Parameterized
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-Representations of Hybrid Automata

We now describe hybrid automata using
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formulas, and define parameterization
and perturbations on them.

A hybrid system is a tuple $H=\langle X$, $Q$, $\mathsf{flow}$,
$\mathsf{guard}$, $\mathsf{reset}$, $\mathsf{inv}$, $\mathsf{init}\rangle$
where $X\subseteq\mathbb{R}^{n}$ specifies the range of the continuous
variables $\vec{x}$ of the system. $Q=\\{q_{0},...,q_{m}\\}$ is a finite set
of discrete control modes. $\mathsf{flow}\subseteq Q\times
X\times\mathbb{R}\times X$ specifies the continuous dynamics for each mode.
The $\mathsf{flow}$ predicate is usually defined either as explicit mappings
from $\vec{a}_{0}$ and $t$ to $\vec{a}_{t}$, or as solutions of systems of
differential equations/inclusions that specify the derivative of $\vec{x}$
over time. $\mathsf{jump}\subseteq Q\times X\times Q\times X$ specifies the
jump conditions between modes. $\mathsf{inv}\subseteq Q\times X$ defines the
invariant conditions for the system to stay in a control mode.
$\mathsf{init}\subseteq Q\times X$ defines the set of initial configurations
of the system. Without loss of generality we always assume that $q_{0}$ is the
only intial mode, and $\mathsf{init}_{q_{0}}\subseteq X$ denotes the initial
values for the continuous variables.

######  Definition 6
($\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-Representations)

Let $H=\langle X$, $Q$, $\mathsf{flow}$, $\mathsf{jump}$, $\mathsf{inv}$,
$\mathsf{init}\rangle$ be an $n$-dimensional hybrid automaton. Let
$\mathcal{F}$ be a set of real functions, and
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$ the corresponding first-order
language. We say that $H$ has an
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-representation, if for every
$q,q^{\prime}\in Q$, there exists quantifier-free
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formulas

$\phi^{q}_{\mathsf{flow}}(\vec{x},\vec{x}_{0},t),\phi^{q\rightarrow
q^{\prime}}_{\mathsf{jump}}(\vec{x},\vec{x}^{\prime}),\phi^{q}_{\mathsf{inv}}(\vec{x}),\phi^{q}_{\mathsf{init}}(\vec{x})$

such that for all $\vec{a},\vec{a}^{\prime}\in\mathbb{R}^{n}$,
$t\in\mathbb{R}$:

  * •

$\mathbb{R}\models\phi^{q}_{\mathsf{flow}}(\vec{a},\vec{a}^{\prime},t)$ iff
$(q,\vec{a},\vec{a}^{\prime},t)\in\mathsf{flow}$.

  * •

$\mathbb{R}\models\phi^{q\rightarrow
q^{\prime}}_{\mathsf{jump}}(\vec{a},\vec{a}^{\prime})$ iff
$(q,q^{\prime},\vec{a},\vec{a}^{\prime})\in\mathsf{jump}$.

  * •

$\mathbb{R}\models\phi^{q}_{\mathsf{inv}}(\vec{a})$ iff
$(q,\vec{a})\in\mathsf{inv}.$

  * •

$\mathbb{R}\models\phi^{q}_{\mathsf{init}}(\vec{a})$ iff $q=q_{0}$ and
$\vec{a}\in\mathsf{init}_{q_{0}}$.

We can write $H=\langle
X,Q,\phi_{\mathsf{flow}},\phi_{\mathsf{jump}},\phi_{\mathsf{inv}},\phi_{\mathsf{init}}\rangle$
to emphasize that $H$ is $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-represented.
But from now on we simply write
$\mathsf{flow},\mathsf{jump},\mathsf{inv},\mathsf{init}$ to denote these logic
formulas, so that we can use $H=\langle
X,Q,\mathsf{flow},\mathsf{jump},\mathsf{inv},\mathsf{init}\rangle$ directly to
denote the $\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-representation of $H$.

######  Definition 7 (Computable Representation)

We say a hybrid automaton $H$ has a computable representation, if $H$ has an
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-representation, where $\mathcal{F}$
is an arbitrary set of computable functions.

Combining continuous and discrete behaviors, the trajectories of hybrid
systems are piecewise continuous. This motivates a two-dimensional structure
of time, with which we can keep track of both the discrete changes and the
duration of each continuous flow.

######  Definition 8 (Hybrid Time Domain)

A hybrid time domain $T$ is a subset of $\mathbb{N}\times\mathbb{R}$ of the
form $T_{m}=\\{(i,t):i<m\mbox{ and }t\in[t_{i},t_{i}^{\prime}]\mbox{ or
}[t_{i},+\infty)\\},$ where $m\in\mathbb{N}\cup\\{+\infty\\}$,
$\\{t_{i}\\}_{i=0}^{m}$ is an increasing sequence in $\mathbb{R}^{+}$,
$t_{0}=0$, and $t_{i}^{\prime}=t_{i+1}$.

We write the set of all hybrid time domains as $\mathbb{H}$.

######  Definition 9 (Hybrid Trajectories)

Suppose $X\subseteq\mathbb{R}^{n}$ and $T_{m}$ is a hybrid time domain. A
hybrid trajectory is any continuous function $\xi:T_{m}\rightarrow X.$

We write $\Xi_{X}$ to denote the set of all possible hybrid trajectories from
$\mathbb{H}$ to $X$. We can now define trajectories of a given hybrid
automaton. The intuition behind the following definition is straightforward.
The labeling function $\sigma_{\xi}^{H}(i)$ is used to map a step $i$ to the
corresponding discrete mode in $H$. In each mode, the system flows
continuously following the dynamics defined by
$\mathsf{flow}(q,\vec{x}_{0},t)$. Note that $(t-t_{k})$ is the actual duration
in the $k$-th mode. When a switch between two modes is performed, it is
required that $\xi(k+1,t_{k+1})$ is updated from the exit value
$\xi(k,t_{k}^{\prime})$ in the previous mode, following the jump conditions.

######  Definition 10 (Trajectories of a Hybrid Automaton)

Let $H$ be a hybrid automaton, and $\xi:T_{m}\rightarrow X$ a hybrid
trajectory. We say that $\xi:T_{m}\rightarrow X$ is a trajectory of $H$ of
discrete depth $m$, if there exists a labeling function
$\sigma^{H}_{\xi}:\mathbb{N}\rightarrow Q$ such that:

  * •

$\sigma^{H}_{\xi}(0)=q_{0}$ and
$\mathbb{R}_{\mathcal{F}}\models\mathsf{init}_{q_{0}}(\xi(0,0))$.

  * •

For any $(i,t)\in T_{m}$,
$\mathbb{R}_{\mathcal{F}}\models\mathsf{inv}_{\sigma^{H}_{\xi}(i)}(\xi(i,t))$.

  * •

When $i=0$,
$\mathbb{R}_{\mathcal{F}}\models\mathsf{flow}_{q_{0}}(\xi(0,0),\xi(0,t),t)$.

  * •

When $i=k+1$, where $0<k+1<m$,

$\displaystyle\mathbb{R}_{\mathcal{F}}$ $\displaystyle\models$
$\displaystyle\mathsf{flow}_{\sigma^{H}_{\xi}(k+1)}(\xi(k+1,t_{k+1}),\xi(k+1,t),$
$\displaystyle(t-t_{k+1}))\mbox{ and }$
$\displaystyle\mathbb{R}_{\mathcal{F}}$ $\displaystyle\models$
$\displaystyle\mathsf{jump}_{(\sigma^{H}(k)\rightarrow\sigma^{H}(k+1))}(\xi(k,t_{k}^{\prime}),\xi(k+1,t_{k+1})).$

We write $\llbracket H\rrbracket$ to denote all possible trajectories of $H$.

######  Definition 11 (Reachability Properties)

Let $H$ be a hybrid automaton and $U\subseteq X\times Q$ be a subset of its
state space. Let $U\subseteq X\times Q$ be a subset of the state space of $H$.
$H$ reaches $U$ if there exists $\xi\in\llbracket H\rrbracket$ such that there
exists $t\in\mathbb{R}$ and $n\in\mathbb{N}$ satisfying

$(\xi(t,n),\sigma_{\xi}^{H}(n))\in U.$

Let $H$ be a hybrid system. Parameter synthesis for reachability properties
asks for a set of parameters such that some mode can be reached.

######  Definition 12 (Parameterized Hybrid Automaton)

We say a hybrid automaton $H$ is parameterized by $\vec{p}$, if We say $H$ is
parameterized by $\vec{p}=(p_{1},...,p_{m})$, if

$H(\vec{p})=\langle
X,Q,\mathsf{flow}(\vec{p}),\mathsf{jump}(\vec{p}),\mathsf{inv}(\vec{p}),\mathsf{init}(\vec{p})\rangle,$

where $\vec{p}$ are among the free variables in the
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-representation of $H$.

######  Definition 13 (Parameter Synthesis for Reachability Properties)

Let $H(\vec{p})$ be a hybrid automaton parameterized by variables
$\vec{p}=(p_{1},...,p_{m})$, and $U\subseteq X\times Q$ a subset of its state
space. Thus, the parameter synthesis problem for reachability asks for an
assignment for $\vec{a}\in\mathbb{R}^{m}$ such that $H(\vec{a})$ reaches $U$.

###  III-C Synthesizing Parameters with $\delta$-Decisions

We now show how to encode parameter synthesis problems for
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-represented hybrid systems using
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formulas. Throughout the following
two definitions, let $H=\langle X$, $Q$, $\mathsf{flow}$, $\mathsf{jump}$,
$\mathsf{init}\rangle$ be an $n$-dimensional
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-represented hybrid system with
$|Q|=m$, and $\mathsf{unsafe}$ an
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formula that encodes a subset
$U\subseteq X\times Q$. Let $k\in\mathbb{N}$ and $M\in\mathbb{R}$ be the
bounds on steps and time respectively. Recall that $q_{0}\in Q$ always denotes
the starting mode.

$\mathsf{Reach}_{H,q^{\prime}}^{k}(\vec{x}_{k}^{t})$ defines the states that
$H$ can reach, if after $k$ steps of discrete changes it is in mode
$q^{\prime}$. From there, if $H$ makes a $\mathsf{jump}$ from mode
$q^{\prime}$ to $q$, then the states have to make a discrete change following
$\mathsf{jump}_{q^{\prime}\rightarrow q}(\vec{x}_{k}^{t},\vec{x}_{k+1})$. As
last, in mode $q^{\prime}$, any state $\vec{x}_{k+1}^{t}$ that $H$ can reach
should satisfy the $\mathsf{flow}$ conditions
$\mathsf{flow}_{q}(\vec{x}_{k+1}^{t},\vec{x}_{k+1},t)$ in mode $q$. Note that
after each discrete jump, a new time variable $t_{k}$ is introduced and
independent from the previous ones.

The $(k,M)$-reachability encoding of $H$ and $U$, $\mathsf{Reach}^{k,M}(H,U)$,
is defined as:

$\displaystyle\exists\vec{a}\exists^{X}\vec{x}_{0}\exists^{X}\vec{x}_{0}^{t}\cdots\exists^{X}\vec{x}_{k}\exists^{X}\vec{x}_{k}^{t}\exists^{[0,M]}t_{0}\cdots\exists^{[0,M]}t_{k}$
$\displaystyle\Big{(}\
\mathsf{init}_{q_{0}}(\vec{x}_{0})\wedge\mathsf{flow}_{q_{0}}(\vec{a},\vec{x}_{0},\vec{x}_{0}^{t},t_{0})$
$\displaystyle\wedge\forall^{[0,t_{0}]}t\forall^{X}\vec{x}\;(\mathsf{flow}_{q_{0}}(\vec{a},\vec{x}_{0},\vec{x},t)\rightarrow\mathsf{inv}_{q_{0}}(\vec{a},\vec{x}))$
$\displaystyle\wedge\bigvee_{i=0}^{k-1}\Big{(}\bigvee_{q,q^{\prime}\in
Q}\Big{(}\mathsf{jump}_{q\rightarrow
q^{\prime}}(\vec{a},\vec{x}_{i}^{t},\vec{x}_{i+1})\wedge\mathsf{flow}_{q^{\prime}}(\vec{a},\vec{x}_{i+1},\vec{x}_{i+1}^{t},t_{i+1})$
$\displaystyle\wedge\forall^{[0,t_{0}]}t\forall^{X}\vec{x}\;(\mathsf{flow}_{q^{\prime}}(\vec{a},\vec{x}_{i+1},\vec{x},t)\rightarrow\mathsf{inv}_{q_{0}}(\vec{a},\vec{x})))\Big{)}\Big{)}$
$\displaystyle\wedge\ \mathsf{unsafe}(\vec{a},\vec{x}_{k}^{t})\Big{)}.$

$H$ reaches $U$ in $k$ steps of discrete jumps with time duration less than
$M$ for each state iff $\mathsf{Reach}^{k,M}(H,U)$ is true.

##  IV Applications

###  IV-A Model Calibration and Falsification

The $\delta$-decision problems can be solved using our dReal tool [52].
Parameter estimation of single-mode ODE models can be encoded as SMT formulas
by BioPSy [53] and solved by dReal, while for multi-mode models we ask a
$k$-step reachability question: Is there a set of parameter values using which
the model reaches the goal region in $k$ steps? The dReach tool [54] can
automatically build such reachability formulas from a multi-mode model and a
goal description, which are then verified by dReal. If $\mathsf{unsat}$ is
returned, the model is unfeasible, which means that the model is unable to
satisfy a desired behavior no matter which parameter values are used. This can
be used to reject model hypotheses. For example, we have showed that the
Fenton-Karma model [55] of cardiac cells is unable to reproduce the “spike-
and-dome” morphology of action potential which has been observed in epicardial
cells [37]. On the other hand, if the model is $\delta$-$\mathsf{sat}$, a
witness (i.e. a set of parameter values) is returned. For example, using the
Bueno-Cherry-Fenton model [56], we have identified critical parameter ranges
that can cause cardiac disorders such as tachycardia and fibrillation [37].

###  IV-B Identification of Therapeutic Strategies

The $\delta$-decision procedures can be used to design optimal therapeutic
strategies. For example, Fig. 3 illustrates a multi-mode model of the TBI-
induced signaling network shown in Fig. 1. The starting point Mode 0
corresponds to live cells under no treatment for 24 h after TBI. Mode 1 is the
“point of no return” which leads to cell death. Modes A-E represent live cells
subjected to particular treatments (e.g. Balcalein-induced ferroptosis
inhibition). The jump conditions are defined by the molecular signature. For
example, starting with Mode 0, if the oxidized CL level exceeds a threshold
$\theta_{1}$, the system jumps to Mode A, which means that we deliver the
apoptosis inhibitor JP4-039. The system then evolves according to the ODEs in
Mode A, until another jump condition–e.g. the level of activated RIP3 above
threshold $\theta_{2}$–is reached, which implies the onset of necroptosis, or
Mode B. We then deliver necroptosis inhibitor necrostatin-1. Suppose the
system next jumps back to Mode 0 and stays. The mode path $0\rightarrow
A\rightarrow B\rightarrow 0$ suggests a successful treatment scheme defined by
a set of jump conditions. Note that we also aim to minimize the number of
drugs used (i.e. path length) to avoid potential side effects. Thus, the
problem of determining which drug to deliver at what time, evolves into a
parameter synthesis problem for hybrid automata and can be tackled by
$\delta$-decision procedures. In a “proof-of-concept” study [38], we have used
this approach to identity personalized therapeutic strategies for prostate
cancer patients.

Figure 3: Multi-mode modeling. Mode 0: live cell, no treatment; Mode A:
apoptosis inhibition; Mode B: necroptosis inhibition; Mode C: Ferroptosis
inhibition; Mode D: Pyroptosis inhibition; Mode E: Parthanatos inhibition Mode
1: dead cell.

###  IV-C Stability Analysis

The robustness of biological systems often refers to the consistency of system
behavior in response to small perturbations. For example, cardiac cells filter
out insignificant stimulations to ensure proper functioning in noisy
environments. Using the $\delta$-decision procedures, we can verify this by
checking if the action potential can be successfully triggered by a small
range of stimulation. An $\mathsf{unsat}$ answer returned by dReach will
guarantee that the model is robust to the corresponding stimulation amplitude
[37].

In addition to “time-bounded” robustness, the $\delta$-decision procedures can
also be used to analyze the infinite-time stability of nonlinear dynamical
systems [57, 58]. Such stability–often referred as structural stability–exists
in the systems that have unique globally asymptotically stable steady states.
A standard way of verifying this property is to find a Lyapunov function [59]
that provides theoretical guarantees on qualitative behavior of the system.
Lyapunov-enable analysis has been applied to mass-action law based kinetic
models of signaling networks such as T-cell kinetic proofreading and ERK
signaling [60]. Our $\delta$-decision procedures enable the Lyapunov stable
analysis for systems with non-polynomial nonlinearity in two ways: (i) Given a
template function, we can synthesize a Lyapunov function by solving
$\exists\forall$-formulas that encode the corresponding non-convex, multi-
objective and disjunctive optimization problem [57]; (ii) We can provide a
sound and relative-complete proof system for induction rules that robustify
the standard notions of Lyapunov functions [58].

##  V Conclusion and Future Works

We have presented a model checking-enabled framework for analyzing systems
biology models with the help of $\delta$-decision procedures. We used the
$\mathcal{L}_{\mathbb{R}_{\mathcal{F}}}$-formulas to describe parameterized
hybrid automata and encode parameter synthesis problems. We employed the
$\delta$-decision procedures to perform bounded model checking and obtain
parameters using which the model can satisfy desired properties. We have
showed that the $\delta$-decision procedures can be used to estimate unknown
model parameters, reject model hypothesis, analyze the model’s Lyapunov
stability, as well as design combination therapies. These tasks can be
assembled as a unified workflow for understanding the mechanism of diseases
and identifying therapeutic options. Our preliminary studies on prostate
cancer, cardiac disorders, and radiation diseases have shown promising
results.

An interesting direction is extending our method for probabilistic settings to
address the inherent variability in biological systems. To cope with the model
complexity, an idea is to approximate the hybrid system as a multi-mode
network of DBNs by extending the approximation technique we have developed for
a single system of ODEs [5]. We plan to explore this in our future work. In
this respect, the Chow-Liu tree representation-based approximation scheme used
in [61] promises to offer helpful pointers.

## Acknowledgment

This work is supported by National Institutes of Health grants P01DK096990,
U19AI068021, and P41GM103712.

## References

  * [1] P. K. Sorger and B. Schoeberl, “An expanding role for cell biologists in drug discovery and pharmacology,” _Mol. Biol. Cell_ , vol. 23, no. 21, pp. 4162–4164, 2012.
  * [2] B. Liu and P. S. Thiagarajan, “Modeling and Analysis of Biopathways Dynamics,” _J. Bioinf. Comput. Biol._ , vol. 10, no. 4, p. 1231001, 2012\.
  * [3] B. Liu, P. Thiagarajan, and D. Hsu, “Probabilistic approximations of signaling pathway dynamics,” in _CMSB’09_ , 2009, pp. 251–265.
  * [4] B. Liu, D. Hsu, and P. Thiagarajan, “Probabilistic approximations of ODEs based bio-pathway dynamics,” _Theor. Comput. Sci._ , vol. 412, no. 21, pp. 2188–2206, 2011.
  * [5] B. Liu _et al._ , “Approximate probabilistic analysis of biopathway dynamics,” _Bioinformatics_ , vol. 28, no. 11, pp. 1508–1516, 2012.
  * [6] A. Hagiescu _et al._ , “GPU code generation for ODE-based applications with phased shared-data access patterns,” _ACM Trans. Archit. Code Optim_ , vol. 10, no. 4, pp. 1–19, 2013.
  * [7] S. K. Palaniappan, S. Akshay, B. Liu, B. Genest, and P. Thiagarajan, “A hybrid factored frontier algorithm for dynamic bayesian networks with a biopathways application,” _IEEE/ACM Trans. Comput. Biol. Bioinform._ , vol. 9, no. 5, pp. 1352–1365, 2012.
  * [8] S. K. Palaniappan, M. Pichené, G. Batt, E. Fabre, and B. Genest, “A look-ahead simulation algorithm for dbn models of biochemical pathways,” in _HSB’16_ , 2016, pp. 3–19.
  * [9] B. Liu _et al._ , “A computational and experimental study of the regulatory mechanisms of the complement system,” _PLoS Comput. Biol._ , vol. 7, no. 1, 2011.
  * [10] S. K. Palaniappan, F. Bertaux, M. Pichené, E. Fabre, G. Batt, and B. Genest, “Abstracting the dynamics of biological pathways using information theory: a case study of apoptosis pathway,” _Bioinformatics_ , vol. 33, no. 13, pp. 1980–1986, 2017.
  * [11] S. K. Palaniappan, B. M. Gyori, B. Liu, D. Hsu, and P. S. Thiagarajan, “Statistical model checking based calibration and analysis of bio-pathway models,” in _CMSB’13_ , vol. 8130 LNCS, 2013, pp. 120–134.
  * [12] B. Liu, B. M. Gyori, and P. Thiagarajan, “Statistical model checking-based analysis of biological networks,” in _Automated Reasoning for Systems Biology and Medicine_ , 2019, pp. 63–92.
  * [13] R. Ramanathan, Y. Zhang, J. Zhou, B. M. Gyori, W.-F. Wong, and P. Thiagarajan, “Parallelized parameter estimation of biological pathway models,” in _HSB’15_ , 2015, pp. 37–57.
  * [14] B. Liu and J. R. Faeder, “Parameter estimation of rule-based models using statistical model checking,” in _BIBM’16_ , 2016, pp. 1453–1459.
  * [15] Q. Wang, N. Miskov-Zivanov, B. Liu, J. R. Faeder, M. Lotze, and E. M. Clarke, “Formal modeling and analysis of pancreatic cancer microenvironment,” in _CMSB’16_ , 2016, pp. 289–305.
  * [16] B. M. Gyori, B. Liu, S. Paul, R. Ramanathan, and P. Thiagarajan, “Approximate probabilistic verification of hybrid systems,” in _HSB’15_ , 2015, pp. 96–116.
  * [17] B. Liu, Q. Liu, S. Palaniappan, I. Bahar, P. S. Thiagarajan, and J. L. Ding, “Innate Immune Memory and Homeostasis May Be Conferred Through TLR3-TLR7 Pathway Crosstalk,” _Sci. Signal._ , vol. 9, no. 436, p. ra70, 2016.
  * [18] B. Liu, D. Bhatt, Z. N. Oltvai, J. S. Greenberger, and I. Bahar, “Significance of p53 dynamics in regulating apoptosis in response to ionizing radiation, and polypharmacological strategies.” _Sci. Rep._ , vol. 4, p. 6245, 2014.
  * [19] B. Liu _et al._ , “Quantitative assessment of cell fate decision between autophagy and apoptosis,” _Sci. Rep._ , vol. 7, no. 1, pp. 1–14, 2017.
  * [20] V. E. Kagan _et al._ , “Oxidized arachidonic and adrenic pes navigate cells to ferroptosis,” _Nat. Chem. Biol._ , vol. 13, no. 1, p. 81, 2017.
  * [21] A. A. Kapralov _et al._ , “Redox lipid reprogramming commands susceptibility of macrophages and microglia to ferroptotic death,” _Nat. Chem. Biol._ , vol. 16, no. 3, pp. 278–290, 2020.
  * [22] J. Steinman _et al._ , “Improved total-body irradiation survival by delivery of two radiation mitigators that target distinct cell death pathways,” _Radiat. Res_ , vol. 189, no. 1, pp. 68–83, 2018.
  * [23] S. Thermozier _et al._ , “Radioresistance of serpinb3a-/- mice and derived hematopoietic and marrow stromal cell lines,” _Radiat. Res._ , vol. 192, no. 3, pp. 267–281, 2019.
  * [24] ——, “Anti-ferroptosis drug enhances total-body irradiation mitigation by drugs that block apoptosis and necroptosis,” _Radiat. Res._ , p. [Epub ahead of print], 2020.
  * [25] R. Goebel, R. G. Sanfelice, and A. R. Teel, “Hybrid dynamical systems,” _IEEE Contr. Syst. Mag._ , vol. 29, no. 2, pp. 28–93, 2009.
  * [26] T. A. Henzinger, “The theory of hybrid automata,” in _LICS’96_ , Washington, DC, USA, 1996, pp. 278–292.
  * [27] R. Ghosh and C. Tomlin, “Symbolic reachable set computation of piecewise affine hybrid automata and its application to biological modelling: Delta-notch protein signalling,” _IET Syst. Biol._ , vol. 1, no. 1, pp. 170–183, 2004.
  * [28] P. Ye, E. Entcheva, S. Smolka, and R. Grosu, “Modelling excitable cells using cycle-linear hybrid automata,” _IET Syst. Biol._ , vol. 2, no. 1, pp. 24–32, 2008.
  * [29] K. Aihara and H. Suzuki, “Theory of hybrid dynamical systems and its applications to biological and medical systems,” _Philos. Trans. R. Soc. London, Ser. A_ , vol. 368, no. 1930, pp. 4893–4914, 2010.
  * [30] M. Antoniotti, B. Mishra, C. Piazza, A. Policriti, and M. Simeoni, “Modeling cellular behavior with hybrid automata: Bisimulation and collapsing,” in _CMSB’03_ , 2003, pp. 57–74.
  * [31] P. Lincoln and A. Tiwari, “Symbolic systems biology: Hybrid modeling and analysis of biological networks,” in _HSCC’04_ , 2004, pp. 660–672.
  * [32] V. Baldazzi, P. T. Monteiro, M. Page, D. Ropers, J. Geiselmann, and H. De Jong, “Qualitative analysis of genetic regulatory networks in bacteria,” in _Understanding the Dynamics of Biological Systems_ , 2011, pp. 111–130.
  * [33] E. Clarke, A. Fehnker, Z. Han, B. Krogh, O. Stursberg, and M. Theobald, “Verification of hybrid systems based on counterexample-guided abstraction refinement,” in _TACAS’03_ , 2003, pp. 192–207.
  * [34] T. A. Henzinger and P. W. Kopke, “Discrete-time control for rectangular hybrid automata,” _Theor. Comput. Sci._ , vol. 221, no. 1, pp. 369–392, 1999.
  * [35] M. Agrawal, F. Stephan, P. Thiagarajan, and S. Yang, “Behavioural approximations for restricted linear differential hybrid automata,” in _HSCC’06_ , 2006, pp. 4–18.
  * [36] S. Gao, J. Avigad, and E. M. Clarke, “Delta-complete decision procedures for satisfiability over the reals,” in _IJCAR’12_ , 2012, pp. 286–300.
  * [37] B. Liu, S. Kong, S. Gao, P. Zuliani, and E. M. Clarke, “Parameter synthesis for cardiac cell hybrid models using $\delta$-decisions,” in _CMSB’14_ , vol. 8859 LNCS, 2014, pp. 99–113.
  * [38] ——, “Towards personalized prostate cancer therapy using delta-reachability analysis,” in _HSCC’15_ , 2015, pp. 227–232.
  * [39] B. Liu, S. Kong, S. Gao, and E. Clarke, “Parameter identification using delta-decisions for biological hybrid systems,” Carnegie Mellon Univ., Tech. Rep. CMU-CS-13-136, 2013.
  * [40] L. Bortolussi and A. Policriti, “Hybrid systems and biology,” in _Formal Methods for Computational Systems Biology_ , 2008, pp. 424–448.
  * [41] R. Alur, “Formal verification of hybrid systems,” in _EMSOFT’ 11, Taipei, Taiwan_ , 2011, pp. 1–6.
  * [42] E. M. Clarke, J. R. Faeder, C. J. Langmead, L. A. Harris, S. K. Jha, and A. Legay, “Statistical model checking in BioLab: Applications to the automated analysis of T-Cell receptor signaling pathway,” in _CMSB’08_ , 2008, pp. 231–250.
  * [43] N. Chabrier-Rivier, M. Chiaverini, V. Danos, F. Fages, and V. SchÃ¤chter, “Modeling and querying biomolecular interaction networks,” _Theor. Comput. Sci._ , vol. 325, no. 1, pp. 25–44, sep 2004.
  * [44] M. Kwiatkowska, G. Norman, and D. Parker, “Using probabilistic model checking in systems biology,” _SIGMETRICS Perform. Eval. Rev._ , vol. 35, no. 4, pp. 14–21, 2008.
  * [45] C. Li, M. Nagasaki, C. H. Koh, and S. Miyano, “Online model checking approach based parameter estimation to a neuronal fate decision simulation model in Caenorhabditis elegans with hybrid functional Petri net with extension,” _Mol. Biosyst._ , vol. 11(Suppl 7), pp. 1–13, 2010.
  * [46] G. Batt, C. belta, and R. Weiss, “Temporal logic analysis of gene networks under parameter uncertainty,” _IEEE T. Automat. Contr._ , vol. 53, pp. 215–229, 2008.
  * [47] A. Donze, G. Clermont, and C. J. Langmead, “Parameter synthesis in nonlinear dynamical systems: Application to systems biology,” _J Comput. Biol._ , vol. 17, no. 3, pp. 325–336, 2010.
  * [48] A. Khalid and S. K. Jha, “Calibration of rule-based stochastic biochemical models using statistical model checking,” in _BIBM’18_. IEEE, 2018, pp. 179–184.
  * [49] A. Pnueli, “The temporal logic of programs,” in _SFCS’77_ , 1977, pp. 46–57.
  * [50] E. M. Clarke, “The birth of model checking,” in _25 Years of Model Checking_ , 2008, pp. 1–26.
  * [51] E. M. Clarke Jr, O. Grumberg, D. Kroening, D. Peled, and H. Veith, _Model Checking_. MIT Press, 2018.
  * [52] S. Gao, S. Kong, and E. M. Clarke, “dReal: An SMT solver for nonlinear theories of reals,” in _CADE’13_ , 2013, pp. 208–214.
  * [53] C. Madsen, F. Shmarov, and P. Zuliani, “Biopsy: an smt-based tool for guaranteed parameter set synthesis of biological models,” in _CMSB’15_ , 2015, pp. 182–194.
  * [54] S. Kong, S. Gao, W. Chen, and E. Clarke, “dreach: $\delta$-reachability analysis for hybrid systems,” in _TACAS’15_ , 2015, pp. 200–205.
  * [55] F. Fenton and A. Karma, “Vortex dynamics in 3D continuous myocardium with fiber rotation: filament instability and fibrillation,” _Chaos_ , vol. 8, pp. 20–47, 1998.
  * [56] A. Bueno-Orovio, E. M. Cherry, and F. H. Fenton, “Minimal model for human ventricular action potentials in tissue,” _J. Theor. Biol._ , vol. 253, pp. 544–560, 2008.
  * [57] S. Kong, A. Solar-Lezama, and S. Gao, “Delta-decision procedures for exists-forall problems over the reals,” in _CAV’18_ , 2018, pp. 219–235.
  * [58] S. Gao _et al._ , “Numerically-robust inductive proof rules for continuous dynamical systems,” in _CAV’19_ , 2019, pp. 137–154.
  * [59] E. D. Sontag, _Mathematical control theory: deterministic finite dimensional systems_ , 2nd ed. Springer, New York, 1998.
  * [60] M. A. Al-Radhawi, D. Angeli, and E. D. Sontag, “A computational framework for a Lyapunov-enabled analysis of biochemical reaction networks,” _PLoS Comput. Biol._ , vol. 16, no. 2, p. e1007681, 2020.
  * [61] M. Pichené, S. K. Palanniappan, E. Fabre, and B. Genest, “Modeling variability in populations of cells using approximated multivariate distributions,” _IEEE/ACM Trans. Comput. Biol. Bioinform._ , pp. 1–12, 2019\.


# Linking the Dynamic PicoProbe Analytical Electron-Optical Beam Line /
Microscope to Supercomputers

Alexander Brace1,2, Rafael Vescovi1, Ryan Chard1, Nickolaus D. Saint2, Arvind
Ramanathan1,2∗, Nestor J. Zaluzec3,4, Ian Foster1,2∗ 1Data Science and
Learning Division, Argonne National Laboratory, Lemont, IL, USA. 2Computer
Science Department, University of Chicago, Chicago, IL, USA. 3Photon Sciences
Division, Argonne National Laboratory, Lemont, IL, USA. 4Physical Sciences
Division, University of Chicago, Chicago, IL, USA.
Contact authors: {ramanathana<EMAIL_ADDRESS>
(2018; 20 February 2007; 12 March 2009; 5 June 2009)

###### Abstract.

The Dynamic PicoProbe at Argonne National Laboratory is undergoing upgrades
that will enable it to produce up to 100s of GB of data per day. While this
data is highly important for both fundamental science and industrial
applications, there is currently limited on-site infrastructure to handle
these high-volume data streams. We address this problem by providing a
software architecture capable of supporting large-scale data transfers to the
neighboring supercomputers at the Argonne Leadership Computing Facility. To
prepare for future scientific workflows, we implement two instructive use
cases for hyperspectral and spatiotemporal datasets, which include: (i) off-
site data transfer, (ii) machine learning/artificial intelligence and
traditional data analysis approaches, and (iii) automatic metadata extraction
and cataloging of experimental results. This infrastructure supports expected
workloads and also provides domain scientists the ability to reinterrogate
data from past experiments to yield additional scientific value and derive new
insights.

automated science, data flow, HPC, AI, ML

††copyright: acmcopyright††journalyear: 2018††doi:
XXXXXXX.XXXXXXX††conference: Make sure to enter the correct conference title
from your rights confirmation emai; June 03–05, 2018; Woodstock, NY††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Applied computing Physical sciences
and engineering

##  1\. Introduction

Experimental facilities around the world rely on computational infrastructure
to support scientific discovery. Such infrastructure is present at all levels
of the “experimental stack”, including: precise control of instrumentation for
atomic-scale measurement, capturing and analyzing big data in real-time, and
publishing such data for the broader research community to access. An
increasingly important design pattern within this paradigm, “closing the
loop,” seeks to tighten the gap between experiment and computation, thereby
increasing the efficiency of the research process and accelerating the rate of
discovery. Machine learning and artificial intelligence (ML/AI) play a lead
role in this pattern as the underlying computational agents that, in many
scenarios, are able to optimize a set of experimental measurements toward an
objective. To achieve this vision, there is a need for robust, open-source,
modular software components to realize end-to-end experimental workflows and
open the door for computationally mediated science.

In this work, we describe our approach to developing such infrastructure for
the Dynamic PicoProbe Analytical Electron-optical Beam Line / Microscope (Sec.
2.1) at Argonne National Laboratory (ANL). The Dynamic PicoProbe is undergoing
a set of upgrades, and upon completion is expected to produce 100s of GB of
data per day during steady-state operation. In the long term, future state-of-
the-art detectors (which will further extend scientific capabilities) will
generate up to 65 GB of data per second ($\approx$200 TB/hour). To prepare for
these intensive data streams, we employed Globus automation services (Chard et
al., 2023; Vescovi et al., 2022) to develop a pair of data flows for
transferring experimental data (hyperspectral and spatiotemporal images) from
the Dynamic PicoProbe site to the Polaris supercomputer at the Argonne
Leadership Computing Facility (ALCF) for analysis and cataloging, as
illustrated in Fig. 1. In addition to transferring data into a permanent store
at ALCF, we provide a simple access portal for researchers to view
experimental analyses to help guide the next set of experimental measurements
and easily share their findings. We leveraged the Globus Search service and
the Django Globus Portal Framework (DGPF) (Saint et al., 2023) to make data
Findable, Accessible, Interoperable, and Reusable (FAIR) (Wilkinson et al.,
2016).

Figure 1. A high-level vision to support computationally mediated science at
the Dynamic PicoProbe. (1) High-dimensional data is generated by the Dynamic
PicoProbe. (2) Data is transferred from on-site computers to more powerful
computing clusters (e.g., ALCF’s Polaris) to perform hyperspectral analysis
using traditional approaches and metadata tracking to enable scientists to
later reinterrogate data. (3) The hyperspectral data is used as input to ML/AI
approaches to: (i) discover interesting patterns in the data, (ii)
characterize the data distribution, (iii) segment and detect features in the
data to assist in calibrating measurement, (iv) perform error correction by
alerting the Dynamic PicoProbe operator to calibration problems. (4) The data
analysis and metadata are synthesized into an actionable summary to assist
domain scientists in performing the measurement(s) (i.e., experiment) of
interest. In this work, we present steps (2) and (3) as different flow use
cases, rather than a single unified flow.

This work presents computational infrastructure to facilitate automated data
transfer from the Dynamic PicoProbe to the Argonne Leadership Computing
Facility (ALCF) in order to:

  * •

Provide a crucial data storage solution for experimental data that would
quickly overwhelm on-site computing resources.

  * •

Use high performance computing to (i) analyze hyperspectral data, and (ii)
leverage ML/AI on spatiotemporal data streams to extract/summarize
scientifically meaningful information from high-volume and high-velocity data
streams.

  * •

Produce a FAIR search index and user portal to catalog experiment metadata and
data products to support scientific campaigns over extended durations and with
multiple users.

Our code, ML/AI model, and datasets are open-source and freely available on
GitHub111https://github.com/ramanathanlab/PicoProbeDataFlow.

##  2\. Methods

###  2.1. Dynamic PicoProbe

Upgrades to the Dynamic PicoProbe (Zaluzec, 2021) at ANL will enable multi-
modal, multi-dimensional, and in-situ characterization of dynamic events at
interfaces in environmental media. Its capabilities will include temporally
resolved (¡ 3 ms) (Scheres, 2012) and sub-atomic ($\sim$ 50 pm)
(Venkatakrishnan et al., 2014) hyperspectral imaging during in-situ/operando
operations in high-vacuum, cryogenic, liquid, and/or gaseous environments of
macromolecular/ionic species of beam-sensitive and soft/hard matter systems.
The Dynamic PicoProbe features: (i) a 30-300 kV monochromated, aberration-
corrected 50 pm electron probe; (ii) sub-atomic imaging capabilities; (iii)
high energy resolution electron spectroscopy (¡ 30 meV); and (iv) the
$X_{PAD}$, the world’s highest collection efficiency hyperspectral x-ray
detector array ($\sim$4.5 sR).

The Dynamic PicoProbe is controlled by a host computer running Windows 10.
Commands are issued through a GUI and haptic control panel. Data from the
instrument are relayed back to the host computer for interactive control and
visualization from four Linux and two Windows 10 systems that control data
acquisition. Windows 10 and macOS user workstations facilitate data
processing, external data transfers, and data backups of user-curated images
and data products (e.g., hyperspectral images). Currently, user machines are
equipped with a 1 Gbps switch that handles external data transfers. Upgrades
are underway to route data directly from the data acquisition system to the
Argonne National Laboratory backbone, which runs at up to 200 Gbps on-site.

###  2.2. Data Flow Infrastructure

This section describes the computational infrastructure that we use to
transfer, analyze, and publish experimental data in near real-time. Each data
flow in this work comprises three distinct processing steps: (i) Data Transfer
with Globus, (ii) Data Analysis, whereby data products are analyzed, plots are
produced, and experiment metadata are extracted, and (iii) Data Publication,
in which the generated plots and experiment metadata are published to a Globus
Search index. We use the term “flow” as a shorthand to describe a data flow in
which multiple stages of computation run serially across heterogeneous
resources and locations. We implement our flows in Python by using the Globus
Architecture for Data-Intensive Experimental Research (Gladier) software
package (Vescovi et al., 2022).

####  2.2.1. Data Transfer

Before a data flow can begin, new data produced from an experiment must be
automatically recognized and used to invoke the flow. While the scientific use
cases highlighted here (Sec. 3.1, 3.2) focus on user-curated data files, we
can apply the same design principles to process data files written directly by
the experimental instrument software. To support automatic data transfers, we
developed a cross-compatible Python application for Windows 10, macOS, and
Linux that uses the watchdog package (Mangalapilly, 2023) to start a new flow
when files are created on the user machine (Sec. 2.1). Our application is very
lightweight as the task logic, orchestration, and fault tolerance are managed
by Gladier/Globus automation services. This software stack allows scaling the
number of concurrent flows (as supported by the available networking
infrastructure) to keep pace with the data-velocity. We also provide an
automatic checkpointing mechanism to avoid undesired flow repeats in cases
where a user needs to resume experimentation after interruption, e.g., if the
user computer needs to be rebooted or the user resumes a set of experiments on
a subsequent day.

When a new file is detected, the Python application starts a Globus flow. Upon
flow start, files are transferred from the user computer to ALCF’s Eagle
storage system, an O(100PB) Lustre file system. In our example use cases, the
files are written in the Electron Microscopy Dataset (EMD) format, a subset of
the Hierarchical Data Format version 5 (HDF5) format that efficiently stores
high dimensional microscopy data (including hyperspectral images and
spatiotemporal images) in a standardized, efficient, binary format. Provisions
are also incorporated to use other cross-platform formats such as the proposed
ISO standard HMSA format (Torpy et al., 2019), as well as parallel formats
used in the scientific community. The data is moved by using the Globus
Transfer service, a cloud-hosted solution for copying data rapidly and
reliably between Globus Connect endpoints. Transfer leverages the OAuth-based
Globus Auth to identify and authenticate users to ensure data is moved
securely.

####  2.2.2. Data Analysis

Once the EMD files arrive on the Eagle file system, two computational steps
are performed: (i) image processing, and (ii) experiment metadata extraction.

For image processing, we employ Globus Compute (Chard et al., 2020), a
federated function-as-a-service platform for secure and reliable remote
computation, to request a compute node on ALCF’s Polaris supercomputer and
thus avoid overwhelming the login nodes. The Globus Compute model employs
user-deployed endpoint agents on remote resources to perform tasks. A user may
submit Python functions for execution by specifying the function body,
arguments, and the endpoint on which the code is to be executed. The Globus
Compute service securely routes the task to the endpoint, where it may either
provision batch resources or perform the task locally before results are
returned to the user via the Compute service. In our case, the endpoint is
configured to acquire compute nodes on the Polaris supercomputer by using the
PBS scheduler.

Next, the EMD file is parsed to extract experiment metadata by using the
HyperSpy Python package (de la Peña et al., 2022). The metadata includes
sample collection date and time; acquisition instrument (i.e., microscope)
details, such as stage and detector positions, beam energy, and magnification;
and other information, such as software versioning. To increase the end-to-end
efficiency of the flow, we combine metadata extraction and the image
processing steps into a single Globus Compute function which avoids reading
the EMD file twice and minimizes flow orchestration overhead.

####  2.2.3. Data Publication

Data publication is achieved by creating and registering the data and
associated metadata (defined by using an extensible schema based on DataCite
(Brase, 2009)) with a Globus Search index. Globus Search is a cloud-hosted
service that builds on ElasticSearch to enable users to create, populate, and
manage indices of searchable metadata. Search provides a fully-featured free-
text search model along with fine-grained security and access control to
facilitate visibility-filtered query results that restrict data
discoverability to authorized users. The publication process is a light-weight
action that transmits the JSON metadata, extracted during the Data Analysis
step (Sec. 2.2.2), to the Search service, and can be performed on a Polaris
login node.

Metadata and results are then visualized by using a Django Globus Portal
Framework (DGPF) data portal. DGPF combines the Django web framework with
Globus to create a customizable data portal based on the Modern Research Data
Portal, a design pattern for providing secure, scalable, and high performance
access to research data. DGPF portals can be used to dynamically display
records in a Globus Search index while leveraging the trusted authentication
systems to render data and results hosted on remote Globus Connect endpoints.
DGPF catalogs have been used to index millions of datasets consisting of many
terabytes of data. In this work, we enable researchers to search their
experimental data and results by the time and date of the associated
experiment.

##  3\. Results

We present results for two scientific use cases, (i) hyperspectral imaging,
and (ii) spatiotemporal imaging, showcasing the generated results from each
application via a user portal display. We also provide a brief vignette that
illustrates the current performance hurdles and suggests areas for
improvement.

###  3.1. Hyperspectral Imaging Data Flow

The hyperspectral image is processed as a 3-dimensional tensor containing
pixel-width, pixel-height, and hyperspectral data. The hyperspectral data
comprises spectroscopic information about both the atomic elemental
composition present in the sample at a given pixel location, as well as the
electron scattering/image data. We generate a plot for users by taking a sum
along the spectroscopy dimension to compute the intensity of the sample at
each pixel, depicted in Fig. 2.A. We also generate a plot of the entire
sample’s spectrum by summing the image over each of the pixel dimensions, as
shown in Fig. 2.B. This spectrum conveys information about the aggregate
atomic composition present throughout the sample. These plots are
automatically rendered in the DGPF data portal along with the extracted
experimental metadata: see Fig. 2.C.

Figure 2. The DGPF interactive user portal allows researchers to quickly
access experimental results and metadata backed by Globus Search. Here we show
(A) a hyperspectral image of a polyamide organic film treated to capture heavy
metals from water (Song et al., 2019), (B) its corresponding spectrum, and (C)
metadata indicating the microscope settings used to collect the data as well
as the atomic composition of the sample.

###  3.2. Spatiotemporal Imaging Data Flow

As a hierarchical file format, EMD files can store many forms of microscopy
data. In this use case, each EMD file stores a multi-dimensional
spatiotemporal image tensor, where the first axis stores the time dimension
and the inner axes store the pixel-width and pixel-height of the image signal.
Here, we investigate a spatiotemporal image with 600 frames showing the motion
of gold nanoparticles on a carbon background. Notably, an additional
hyperspectral dimension (Sec. 3.1) could be added which would result in a
4-dimensional tensor, vastly increasing the data volume of each file—we leave
this use case to future work.

Figure 3. Spatiotemporal data is preprocessed so that each time step is input
to a fine-tuned YOLOv8 model to track locations of gold nanoparticles. The
bounding boxes (orange) have a confidence score and can be used to count the
number of nanoparticles likely to be in a sample, helping to characterize
changes in the sample as a function of time.

In order to prepare for such data streams, we demonstrate ML/AI approaches to
automatically track scientifically meaningful information about sample
contents. Specifically, we train a YOLOv8 (Jocher et al., 2023) model to
detect and track gold nanoparticles/nanostructures as they move. Before
training a YOLOv8 model, hand-labeled bounding boxes must be drawn around the
prediction targets (in this case, a single label type representing the gold
nanoparticles). As the dataset presented in this work has 600 time steps, we
select every 50th step for hand-labeling with Roboflow (Dwyer et al., 2022),
yielding a total of nine training, three validation, and one testing
640$\times$640 image(s). We further augmented the training set by using
horizontal and vertical flips, as well as random cropping up to 20% maximum
zoom. We then fine-tuned the YOLOv8s model (11.2M parameters) in Google Colab
for 100 epochs with stochastic gradient descent, with a batch size of 16 and a
learning rate of 0.01, on a Tesla T4 GPU. Our model achieves a mean Average
Precision with an Intersection over Union (IoU) range of 50–95% (mAP50-95) of
0.791 on the training set and 0.801 on the validation set.

We employ the fine-tuned YOLOv8 model for efficient inference within the
spatiotemporal imaging data flow by first converting incoming EMD files to MP4
video format, followed by calling the inference routine in a subprocess. We
performed inference on a Polaris compute node with an NVIDIA A100 GPU and
output an annotated MP4 file containing the predicted gold nanoparticle
bounding boxes, as illustrated in Fig. 3.

###  3.3. Performance Evaluation

To provide a controlled environment for testing our data flow infrastructure,
we employ an application that periodically copies a file into the transfer
directory of the Dynamic PicoProbe user computer to simulate data generation
during an actual experiment. We configure the experiments based on the
approximate time it takes each transfer to complete. Thus, over the course of
an hour, we automatically start a new flow every 30 and 120 seconds for the
hyperspectral and spatiotemporal use cases, respectively. The Globus services
allow parallel flow execution that enables us to start new flows even when
previous ones are still running. We summarize the performance metrics for our
use cases in Table LABEL:tab:performance. Note that the file size in the
hyperspectral use case (91 MB) is much smaller than the spatiotemporal
counterpart (1200 MB), which leads to many more hyperspectral flow runs
completing within the allotted hour. The maximum runtimes are associated with
the first flows, as they have to request a compute node on Polaris and cache
the Python libraries required for analysis. Subsequent flows are able to reuse
nodes already provisioned to the previous flows.

Table 1. Performance measured during independent 1-hour long experiments in which files were transferred from the Dynamic PicoProbe user computer to the Eagle filesystem. Metric | Hyperspectral | Spatiotemporal
---|---|---
Start period (s) | 30 | 120
Transfer volume (MB) | 91 | 1200
Total data transfer (GB) | 6.42 | 21.72
Min flow runtime (s) | 29 | 195
Mean flow runtime (s) | 47 | 224
Max flow runtime (s) | 181 | 274
Median overhead (s) | 19.5 | 45.2
Median overhead (%) | 49.2 | 21.1
Total flow runs | 72 | 18

Figure 4. The itemized runtime statistics (in seconds) for the hyperspectral
flow (A) and the spatiotemporal flow (B) measured over independent 1-hour long
experiments. We use “Active” to denote the time spent actively processing
either the Transfer, Analysis, or Publication steps. The overhead measures the
remainder of the total flow runtime not spent actively processing the steps.

In addition to aggregate flow statistics, we characterize the performance of
the individual flow component steps (Sec. 2.2) to profile and understand
performance bottlenecks. We illustrate the runtime statistics (in seconds) of
the hyperspectral flow in Fig. 4.A, and spatiotemporal flow in Fig. 4.B, as
well as the time spent actively processing tasks versus the overhead. The flow
orchestration overhead is significant at 49.2% of the total median runtime for
the hyperspectral flow and 21.1% for the spatiotemporal flow. This observation
is attributed to an exponential polling backoff policy that starts at 1 second
and doubles up to 10 minutes, which we are working to improve. Outside of
overhead, the file transfer time, which dominates active flow runtime, is
primarily a function of the size of the transferred files and demonstrates the
need to update on-site data transfer capabilities (currently facilitated by a
1 Gbps switch (Sec. 2.1)) to support future detectors, which will produce up
to 65 GB per second. The spatiotemporal compute phase can also be optimized
since the majority of time is spent on converting raw EMD files to MP4 format,
which involves a slow data type casting operation from fp64 to uint8. More
efficient integration with the YOLOv8 algorithm would lead to a substantial
improvement in time-to-solution for spatiotemporal data stream analysis.

##  4\. Related Work

We distinguish in this work between computationally mediated science and
automated science. In the former, automation and intelligent agents augment
the abilities of a human scientist by performing routine tasks, automatic
calibrations, and sharing in a collaborative synergy to accelerate discovery
(Zaluzec, 2003). On the other hand, automated science seeks to create self-
driving laboratories that integrate a variety of instruments under the
direction of an ML/AI agent, which automatically steers an experimental
campaign towards discoveries using minimal human intervention (Vescovi et al.,
2023; King et al., 2009). These directions share a common need for modular
software infrastructure (Balouek-Thomert et al., 2019; Salim et al., 2019;
Ananthakrishnan et al., 2020; Vescovi et al., 2022; Zhang et al., 2021) to
link experimental facilities to high performance computing resources (Wang et
al., 2001; Veseli et al., 2018; Basu et al., 2019; Liu et al., 2021).

ML/AI has been successfully applied to a variety of electron microscopy tasks
including: detecting microstructures (Zhong et al., 2022), image registration
(Luo et al., 2021), pixel-level nanoparticle segmentation (Larsen et al.,
2022), frame-level defect tracking and detection (Shen et al., 2021), and many
other applications such as particle picking (Wagner et al., 2019), automated
labeling (Weber et al., 2018), denoising (Bepler et al., 2020), and super-
resolution reconstruction (Suveer et al., 2019). Treder et al. (Treder et al.,
2022) provide a comprehensive review of deep learning in electron microscopy.
Our work extends these efforts by establishing the infrastructure needed to
bridge microscopes with HPC infrastructure for online integration of
computationally expensive ML/AI and traditional analysis techniques.

##  5\. Conclusion

We have presented software infrastructure that links the Dynamic PicoProbe to
supercomputing resources to (1) provide petabyte-scale data storage, (2)
enable near real-time data analysis by using ML/AI techniques, and (3) build
interactive FAIR data portals for researchers to view results and inform
future experiments. We leverage Globus automation services to implement two
prototypical science use cases (hyperspectral imaging and spatiotemporal
imaging), providing configurable software that decouples data analysis steps
from the limitations of on-site computers. We found transfer times between on-
site data staging systems and the supercomputing facility to be the overall
bottleneck in our data flows. We expect this issue to compound as data tensor
dimensions are introduced to measure additional physical parameters such as
temperature and pressure. As such, active directions for future research
include: (1) on-site hardware upgrades, (2) data compression algorithms, and
(3) optimization of cross-site transfer settings.

###### Acknowledgements.

The Analytical Picoprobe at Argonne National Laboratory (ANL) was developed as
part of CRADA #01300710 between ANL and ThermoFisher Scientific Instruments.
This work is supported by the U. S. Department of Energy, Office of Science,
Office of Basic Energy Sciences, under Contract No. DE-AC02-06CH11357, as well
as the National Science Foundation Major Research Instrumentation (MRI)
Program (NSF DMR-2117896) at the University of Chicago. This work was
supported in part by NSF grants OAC-1835890 and OAC-2004894. We are grateful
to the Globus team for their support. We thank Carla M. Mann for helping to
proofread this manuscript.

## References

  * (1)
  * Ananthakrishnan et al. (2020) Rachana Ananthakrishnan, Kyle Chard, Mike D’Arcy, Ian Foster, Carl Kesselman, Brendan McCollam, Jim Pruyne, Philippe Rocca-Serra, Robert Schuler, and Rick Wagner. 2020.  An open ecosystem for pervasive use of persistent identifiers.  In _Practice and Experience in Advanced Research Computing_. 99–105.
  * Balouek-Thomert et al. (2019) Daniel Balouek-Thomert, Eduard Gibert Renart, Ali Reza Zamani, Anthony Simonet, and Manish Parashar. 2019\.  Towards a computing continuum: Enabling edge-to-cloud integration for data-driven workflows.  _International Journal of High Performance Computing Applications_ 33, 6 (2019), 1159–1174.
  * Basu et al. (2019) Shibom Basu, Jakub W Kaminski, Ezequiel Panepucci, C-Y Huang, Rangana Warshamanage, Meitian Wang, and Justyna Aleksandra Wojdyla. 2019\.  Automated data collection and real-time data analysis suite for serial synchrotron crystallography.  _Journal of Synchrotron Radiation_ 26, 1 (2019), 244–252.
  * Bepler et al. (2020) Tristan Bepler, Kotaro Kelley, Alex J Noble, and Bonnie Berger. 2020.  Topaz-Denoise: General deep denoising models for cryoEM and cryoET.  _Nature Communications_ 11, 1 (2020), 5208.
  * Brase (2009) Jan Brase. 2009\.  DataCite–A global registration agency for research data. In _4th International Conference on Cooperation and Promotion of Information Resources in Science and Technology_. IEEE, 257–261.
  * Chard et al. (2020) Ryan Chard, Yadu Babuji, Zhuozhao Li, Tyler Skluzacek, Anna Woodard, Ben Blaiszik, Ian Foster, and Kyle Chard. 2020\.  FuncX: A federated function serving fabric for science. In _29th International Symposium on High-performance Parallel and Distributed Computing_. 65–76.
  * Chard et al. (2023) Ryan Chard, Jim Pruyne, Kurt McKee, Josh Bryan, Brigitte Raumann, Rachana Ananthakrishnan, Kyle Chard, and Ian T Foster. 2023.  Globus automation services: Research process automation across the space–time continuum.  _Future Generation Computer Systems_ 142 (2023), 393–409.
  * de la Peña et al. (2022) Francisco de la Peña, Eric Prestat, Vidar Tonaas Fauske, Pierre Burdet, Jonas Lähnemann, Petras Jokubauskas, Tom Furnival, Magnus Nord, Tomas Ostasevicius, Katherine E. MacArthur, Duncan N. Johnstone, Mike Sarahan, Joshua Taillon, Thomas Aarholt, pquinn dls, Vadim Migunov, Alberto Eljarrat, Jan Caron, Carter Francis, T. Nemoto, Timothy Poon, Stefano Mazzucco, actions user, Nicolas Tappy, Niels Cautaerts, Suhas Somnath, Tom Slater, Michael Walls, Florian Winkler, and Håkon Wiik Ånes. 2022.  _hyperspy/hyperspy: Release v1.7.3_.  https://doi.org/10.5281/zenodo.7263263
  * Dwyer et al. (2022) B Dwyer, J Nelson, J Solawetz, et al. 2022\.  Roboflow (version 1.0)[software].
  * Jocher et al. (2023) Glenn Jocher, Ayush Chaurasia, and Jing Qiu. 2023\.  _Ultralytics YOLOv8_.  https://github.com/ultralytics/ultralytics
  * King et al. (2009) Ross D King, Jem Rowland, Stephen G Oliver, Michael Young, Wayne Aubrey, Emma Byrne, Maria Liakata, Magdalena Markham, Pinar Pir, Larisa N Soldatova, Andrew Sparkes, Kenneth E. Whelan, and Amanda Clare. 2009.  The automation of science.  _Science_ 324, 5923 (2009), 85–89.
  * Larsen et al. (2022) Matthew Helmi Leth Larsen, Cuauhtémoc Nuñez Valencia, William Bang Lomholdt, Daniel Kelly, Pei Liu, Jakob Birkedal Wagner, Ole Winther, Thomas Willum Hansen, and Jakob Schiøtz. 2022.  Large-scale Automated Analysis of High-Resolution Transmission Electron Microscopy Data Assisted by Deep Learning Neural Networks.  _Microscopy and Microanalysis_ 28, S1 (2022), 2984–2986.
  * Liu et al. (2021) Zhengchun Liu, Ahsan Ali, Peter Kenesei, Antonino Miceli, Hemant Sharma, Nicholas Schwarz, Dennis Trujillo, Hyunseung Yoo, Ryan Coffee, Naoufal Layad, Jana Thayer, Ryan Herbst, Chunhong Yoon, and Ian Foster. 2021\.  Bridging data center AI systems with edge computing for actionable information retrieval. In _3rd Annual Workshop on Extreme-scale Experiment-in-the-Loop Computing_. IEEE, 15–23.
  * Luo et al. (2021) Yanqi Luo, Nestor Zaluzec, Mathew Cherukara, Xiaolan Wu, and Si Chen. 2021.  Real-time image registration via a deep leaning approach for correlative X-ray and electron microscopy.  _Microscopy and Microanalysis_ 27, S1 (2021), 302–304.
  * Mangalapilly (2023) Yesudeep Mangalapilly. 2023\.  watchdog.  https://github.com/gorakhargosh/watchdog.
  * Saint et al. (2023) Nickolaus Saint, Ryan Chard, Rafael Vescovi, Jim Pruyne, Ben Blaiszik, Rachana Ananthakrishnan, Michael E Papka, Kyle Chard, and Ian Foster. 2023.  Active Research Data Management with the Django Globus Portal Framework.  (2023).
  * Salim et al. (2019) Michael Salim, Thomas Uram, J Taylor Childers, Venkatram Vishwanath, and Michael Papka. 2019\.  Balsam: Near real-time experimental data analysis on supercomputers. In _1st Annual Workshop on Large-scale Experiment-in-the-Loop Computing_. IEEE, 26–31.
  * Scheres (2012) Sjors HW Scheres. 2012\.  A Bayesian view on cryo-EM structure determination.  _Journal of Molecular Biology_ 415, 2 (2012), 406–418.
  * Shen et al. (2021) Mingren Shen, Guanzhao Li, Dongxia Wu, Yudai Yaguchi, Jack C Haley, Kevin G Field, and Dane Morgan. 2021.  A deep learning based automatic defect analysis framework for In-situ TEM ion irradiations.  _Computational Materials Science_ 197 (2021), 110560.
  * Song et al. (2019) Xiaohui Song, John W Smith, Juyeong Kim, Nestor J Zaluzec, Wenxiang Chen, Hyosung An, Jordan M Dennison, David G Cahill, Matthew A Kulzick, and Qian Chen. 2019\.  Unraveling the morphology-function relationships of polyamide membranes using quantitative electron tomography.  _ACS Applied Materials & Interfaces_ 11, 8 (2019), 8517–8526.
  * Suveer et al. (2019) Amit Suveer, Anindya Gupta, Gustaf Kylberg, and Ida-Maria Sintorn. 2019.  Super-resolution reconstruction of transmission electron microscopy images using deep learning. In _IEEE 16th International Symposium on Biomedical Imaging_. IEEE, 548–551.
  * Torpy et al. (2019) Aaron Torpy, Mike Kundmann, Nicholas C. Wilson, Colin M. MacRae, and Nestor J. Zaluzec. 2019\.  _HMSA File Format Specification - Version 1.02 (November 2019)_.  Technical Report. Microscopy Society of America.
  * Treder et al. (2022) Kevin P Treder, Chen Huang, Judy S Kim, and Angus I Kirkland. 2022.  Applications of deep learning in electron microscopy.  _Microscopy_ 71, Supplement_1 (2022), i100–i115.
  * Venkatakrishnan et al. (2014) SV Venkatakrishnan, Lawrence F Drummy, Michael Jackson, Marc De Graef, Jeff Simmons, and Charles A Bouman. 2014.  Model-based iterative reconstruction for bright-field electron tomography.  _IEEE Transactions on Computational Imaging_ 1, 1 (2014), 1–15.
  * Vescovi et al. (2022) Rafael Vescovi, Ryan Chard, Nickolaus D Saint, Ben Blaiszik, Jim Pruyne, Tekin Bicer, Alex Lavens, Zhengchun Liu, Michael E Papka, Suresh Narayanan, et al. 2022\.  Linking scientific instruments and computation: Patterns, technologies, and experiences.  _Patterns_ 3, 10 (2022).
  * Vescovi et al. (2023) Rafael Vescovi, Tobias Ginsburg, Kyle Hippe, Doga Ozgulbas, Casey Stone, Abraham Stroka, Rory Butler, Ben Blaiszik, Tom Brettin, Kyle Chard, et al. 2023\.  Towards a Modular Architecture for Science Factories.  _arXiv preprint arXiv:2308.09793_ (2023).
  * Veseli et al. (2018) Siniša Veseli, Nicholas Schwarz, and Collin Schmitz. 2018\.  APS data management system.  _Journal of Synchrotron Radiation_ 25, 5 (2018), 1574–1580.
  * Wagner et al. (2019) Thorsten Wagner, Felipe Merino, Markus Stabrin, Toshio Moriya, Claudia Antoni, Amir Apelbaum, Philine Hagel, Oleg Sitsel, Tobias Raisch, Daniel Prumbaum, Dennis Quentin, Daniel Roderer, Sebastian Tacke, Birte Siebolds, Evelyn Schubert, Tanvir R. Shaikh, Pascal Lill, Christos Gatsogiannis, and Stefan Raunser. 2019.  SPHIRE-crYOLO is a fast and accurate fully automated particle picker for cryo-EM.  _Communications Biology_ 2, 1 (2019), 218.
  * Wang et al. (2001) Yuxin Wang, Francesco De Carlo, Derrick C Mancini, Ian McNulty, Brian Tieman, John Bresnahan, Ian Foster, Joseph Insley, Peter Lane, Gregor von Laszewski, Carl Kesselman, Mei-Hui Su, and Marcus Thiebaux. 2001.  A high-throughput x-ray microtomography system at the Advanced Photon Source.  _Review of Scientific Instruments_ 72, 4 (2001), 2062–2068.
  * Weber et al. (2018) Gunther H Weber, Colin Ophus, and Lavanya Ramakrishnan. 2018\.  _Automated labeling of electron microscopy images using deep learning_.  IEEE.
  * Wilkinson et al. (2016) Mark D Wilkinson, Michel Dumontier, IJsbrand Jan Aalbersberg, Gabrielle Appleton, Myles Axton, Arie Baak, Niklas Blomberg, Jan-Willem Boiten, Luiz Bonino da Silva Santos, Philip E Bourne, et al. 2016\.  The FAIR Guiding Principles for scientific data management and stewardship.  _Scientific Data_ 3, 1 (2016), 1–9.
  * Zaluzec (2021) Nestor Zaluzec. 2021\.  First light on the Argonne PicoProbe and the X-ray perimeter array detector (XPAD).  _Microscopy and Microanalysis_ 27, S1 (2021), 2070–2074.
  * Zaluzec (2003) Nestor J Zaluzec. 2003\.  Computationally mediated experimental science.  _Microscopy and Microanalysis_ 9, S02 (2003), 150–151.
  * Zhang et al. (2021) Qingteng Zhang, Eric M Dufresne, Yasukazu Nakaye, Pete R Jemian, Takuto Sakumura, Yasutaka Sakuma, Joseph D Ferrara, Piotr Maj, Asra Hassan, Divya Bahadur, Subramanian Ramakrishnan, Faisal Khan, Sinisa Veseli, Alec R Sandy, Nicholas Schwarz, and Suresh Narayanan. 2021.  20 $\mu$s-resolved high-throughput X-ray photon correlation spectroscopy on a 500k pixel detector enabled by data-management workflow.  _Journal of Synchrotron Radiation_ 28, 1 (2021), 259–265.
  * Zhong et al. (2022) Xiaoting Zhong, Nestor J Zaluzec, Yu Lin, and Jiadong Gong. 2022\.  Machine Learning Enabled Reproducible Data Analysis for Electron Microscopy.  _Microscopy and Microanalysis_ 28, S1 (2022), 3138–3140.