Experimental Biological Protocols with Formal
Semantics
Alessandro Abate2 , Luca Cardelli1,2 , Marta Kwiatkowska2 , Luca Laurenti2 ,
and Boyan Yordanov1

arXiv:1710.08016v1 [] 22 Oct 2017

1

2

Microsoft Research Cambridge
Department of Computer Science, University of Oxford

Abstract. Both experimental and computational biology is becoming
increasingly automated. Laboratory experiments are now performed automatically on high-throughput machinery, while computational models
are synthesized or inferred automatically from data. However, integration between automated tasks in the process of biological discovery is still
lacking, largely due to incompatible or missing formal representations.
While theories are expressed formally as computational models, existing
languages for encoding and automating experimental protocols often lack
formal semantics. This makes it challenging to extract novel understanding by identifying when theory and experimental evidence disagree due
to errors in the models or the protocols used to validate them. To address
this, we formalize the semantics of a core protocol language as a stochastic hybrid process, which provides a unified description for the models of
biochemical systems being experimented on, together with the discrete
events representing the liquid-handling steps of biological protocols. Such
a representation captures uncertainties in equipment tolerances, making
it a suitable tool for both experimental and computational biologists. We
illustrate how the proposed protocol language can be used for automated
verification and synthesis of laboratory experiments on case studies from
the fields of chemistry and molecular programming.

1

Introduction

The classical cycle of observation, hypothesis formulation, experimentation, and
falsification, which has driven scientific and technical progress since the scientific
revolution, is lately becoming automated in all its separate components. Data
gathering is conducted by high-throughput machinery. Models are automatically
synthesized, at least in part, from data [6,11,4]. Experiments are selected to maximize knowledge acquisition. Laboratory protocols are run under reproducible
and auditable software control. However, integration between these automated
components is lacking. Theories are not placed in the same formal context as the
(coded) protocols that are supposed to test them. Theories talk about changing
in physical quantities, while protocols talk about steps carried out by machines:
neither knows about the other, although they both try to describe the same
process. The consequence is that often it is hard to tell what happened when
experiments and models do not match: was it an error in the model, or an error
in the protocol? Often both the model and the protocol have unknown parameters: do we use the experimental data to fit the model or to fit the protocol?

When most activities are automated, we need a way to answer those questions
that is equally automated.
In this paper, we present a core language to model experimental biological
protocols that gives an integrated description of the protocol and of the underlying molecular process. A basic example of experimental biological protocol is
shown in Figure 1.
From this integrated representation both the model of a phenomenon
(for possibly automated mathematical analysis), and the steps carried out
to test it (for automated execution by
lab equipment) can be separately extracted. This is essential to perform
automated model synthesis and falsification by taking into account uncertainties in the both model structure
and in equipment tolerances. We map
our language into a Piecewise Deterministic Markov Process (PDMP).
That is, a class of Markov stochastic
hybrid processes where the continuous
variables evolve according to ordinary
differential equations (ODEs) and the Fig. 1: Graphical representation of an
discrete variables evolve by means of acid-base titration protocol. The protorandom jumps [12]. The discrete dy- col is initialized with samples A (connamics are used to map the discrete taining H + and Cl− ) and B (containoperations of a lab protocol, while ing N a+ and OH − ). Some fraction of
continuous dynamics model the evo- each sample (p1 and p2 ) is mixed tolution of the physical variables. In our gether and the resulting sample is let to
language physical variables are de- equilibrate for t seconds.
scribed with Chemical Reaction Networks (CRNs), a widely used formalism to model molecular interactions [16].
Our goal is to define a simple core language and focusing on formalizing its
semantics. We then show how our language can easily be extended to collect
observations of the process and to model complicate protocols. Giving to the
language a formal semantics in terms of a PDMP allows us to include in our semantics the uncertainties intrinsic in the discrete operations of an experimental
protocol. Uncertainties for experimental protocols have also been standardized
(standards ISO 17025 and 8655). On examples from chemistry and molecular
programming, we demonstrate how our integrated representation allows one to
perform analysis and synthesis of both the discrete steps of the protocol and of
the underlying physical system.
Related Work Several factors contribute to the growing need for a formalization of experimental protocols in biology. First, better record-keeping of experimental operations is recognized as a step towards tackling the reproducibility

crisis in biology [17]. Second, the emergence of cloud labs [18] creates a need for
precise, machine-readable descriptions of the experimental steps to be executed.
To address these needs, frameworks allowing protocols to be recorded, shared,
and reproduced locally or in a remote lab have been proposed. These frameworks
introduce different programming languages for experimental protocols including
BioCoder [3], Autoprotocol, and Antha [24]. These languages provide expressive, high-level protocol descriptions but consider each experimental sample as
a labelled black-box. This makes it challenging to study a protocol together
with the biochemical systems it manipulates in a common framework. In contrast, we consider a simpler set of protocol operations but capture the details of
experimental samples, enabling us to track properties of chemical species (e.g.
amounts, concentrations, etc) as they react during the execution of a protocol.
This allows us to formalize and verify requirements for the correct execution of
a protocol or to optimize various protocol or system parameters to satisfy these
specifications.

2

Background

We first introduce the formalism of PDMP, which we use to model experimental
protocols. Then, we introduce Chemical Reaction Networks (CRN), which are
used to model the underlying physical process.
2.1 Piecewise Deterministic Markov Process
The syntax of a PDMP is given as follows.
Definition 1 A Piecewise Deterministic Markov Process (PDMP) H is a tuple
H = (Q, d, G, F, Λ, R), where
– Q = {q1 , ..., q|Q| } is the set of discrete modes
– d : Q → N is a map such that Rd(q) , is the state space of the continuous
dynamics for state q. The hybrid state space is defined as D = ∪q∈Q {q} ×
Rd(q)
– G : Q × Rd(q) → {0, 1} is a set of guards
– F : Q × Rd(q) → Rd(q) is a family of vector fields
– Λ : S × Q → R≥0 is an intensity function,
where for (qi , x) ∈ S, qj ∈ Q, we
P
define Λ((qi , x), qj ) = λi,j (x) and qj 6=qi λi,j (x) = λqi (x)
– R : B(S) × S → [0.1] is the reset function, which assigns to each (q, x) ∈ S
a measure R(·, q, x) on (S, B(S)).
Let B(S) be the smallest σ−algebra on S containing all the sets of the form
∪q∈Q {q}×Aq , where Aq is a Borel subeset of Rd(q) . For t ∈ R≥0 , q ∈ Q, x ∈ Rd(q) ,
we call Φ(q, t, x), the solution of the following differential equation
dΦ(q, t, x)
= F (q, Φ(q, t, x)), Φ(q, 0, x) = x.
dt
The solution of a PDMP is a stochastic process Y = (α, X), whose semantics is
classically defined according to the notion of execution (see Definition 2 below)
[13]. In order to introduce such a notion, we define the exit time t∗ (q, x, G) as
t∗ (q, x, G) = inf{t ∈ R≥0 | G(q, Φ(q, t, x)) = 1}

(1)

(

Rt

e− 0 λq (Φ(q,τ,x))dτ
if t < t∗ (q, x, G)
.
and the survival function f (q, t, x) =
0
otherwise.
Here t∗ (q, x, G) represents the first time instant, starting from state (q, x), when
the guard set is reached by a solution of the process; further f (q, t, x) denotes
the probability that the system remains within q, starting from x, at time t
[12], which depends on random arrivals induced by the intensity function Λ. The
semantics of a PDMP is provided next.
Definition 2 (Execution of PDMP H)
Set t := 0
Set (α(0), X(0)) := q0 , x0
While t < ∞
Extract R≥0 -valued random variable T such that
P rob(T > t̄) = f (α(t), t̄, X(t))
∀τ ∈ [t, t + T ) Set (α(τ ), X(τ )) := (α(t), Φ(α(t), τ − t, X(t)))
If t + T < ∞
Extract (α(t + T ), X(t + T )) according to
R(·, (α(t), Φ(α(t), T, X(t)))
End If
Set t := t + T
End While
Let T < t∗ (qi , x̄, G) be the dwelling time in state qi ∈ Q, with x̄ such that
Φ(qi , T, x̄) = x, (qj , A) ∈ B(S). Assume that x is such that G(x, qi ) = 0. Then,
the reset has the following form, which results from a transition that is not due
to the crossing of a guard:
R((qj , A), (qi , x)) = Rc (A|qj , (qi , x)) P

ΛTqi ,qj (Φ, x)

qk 6=qi

where ΛTqi ,qj (Φ, x) =

RT
0

ΛTqi ,qk (Φ, x)

,

(2)

λi,j (Φ(s, qi , x))ds, and

Rc (A|qj , (qi , x)) = P rob(X(T ) ∈ A|α(T ) = qj , (α(0) = qi , X(0) = x))
is the conditional reset of the continuous dynamics.
2.2 Chemical Reaction Networks
A CRN C = (A, R) is a pair of finite sets, where A denotes a set of chemical
species, |A| is its cardinality, and R denotes a set of reactions. A reaction τ ∈ R
is a triple τ = (rτ , pτ , kτ ), where rτ ∈ N|Λ| is the source complex, pτ ∈ N|Λ| is
the product complex and kτ ∈ R>0 is the coefficient associated with the rate of
the reaction.The quantities rτ and pτ represent the stoichiometry of reactants
and products. Given a reaction τ1 = ([1, 0, 1], [0, 2, 0], k1 ) we often refer to it
visually as τ1 : λ1 + λ3 →k1 2λ2 . The net change associated to τ is defined by
υτ = pτ − rτ .

Many models have been introduced to study CRNs [9,7,15,8]. Here we consider the reaction rate equations [15], which describe the time evolution of the
concentration of the species in C, in a sample of temperature T and volume V ,
as follows:
X
dΦ(t)
= F (t) =
υτ · γS (Φ(t), kτ , V, T ),
(3)
dt
τ ∈R

where γS (Φ(t), kτ , V, T )) is the propensity rate, and in case of mass action kinetics we have
Y r
γS (Φ(t), kτ , V, T )) = kτ (T )
ΦSS,τ (t),
S∈Λ

where ΦS and rS,τ are the components of vectors Φ and rτ relative to species S,
and where in kτ (T ) we make explicit the dependence from temperature T .
Definition 3 (Chemical Reaction System) A chemical reaction system (CRS)
C = (A, R, x0 ) is defined as a tuple, where (A, R) is a CRN and x0 ∈ N|Λ|
represents its initial condition.
Example 1. Consider the CRS C = (A, R, x0 ), evolving in a volume V and at
temperature T , where A = {H2 O, N a+ , OH − , Cl− , H + } and R is composed of
the following reactions:
N a+ + OH − + H + + Cl− →k H2 O + N a+ + Cl−
where k = 2.81e−10 is the rates at temperature T = 298 Kelvin. Then, According
to Equation (3), we have that the state of H + is given by the solution of the
following ordinary differential equation:
dH + (t)
= − kN a+ (t)OH − (t)H + (t)Cl− (t),
dt
with H + (0) =

3

x0,H +
V

, where x0,H + is the component of x0 relative to H + .

A Language for Experimental Biological Protocols

In this section we introduce the syntax of a language for modelling experimental
protocols. A formal semantics of the language, based on denotational semantics
[25], is then discussed. We model the physical process underlying a biological
experimental protocol as a CRS. As a consequence, in order to introduce formal
semantics for experimental protocols, we first need to define semantics for a CRS,
which has been only introduced informally in the previous section.
Let S = (R|A| × R≥0 × R≥0 ) be a sample. We define the semantics for a CRS
as follows.
|A|

Definition 4 (CRS Semantics) Let C = (A, R) be a CRN, x0 ∈ R≥0 , V, T ∈
R≥0 be the initial concentration (moles), volume (liters) and temperature (degrees
Kelvin). Call F (V, T ) : R|A| → R|A| the drift at volume V and temperature T

for C. Then, the semantics of the CRS (A, R, x0 ) at volume V , temperature T
and time t, for a time horizon H ∈ R≥0 ∪ {∞},
[[·]] : (CRS × R≥0 × R≥0 ) → R≥0 ∪ {∞} → R≥0 → S
is defined as
[[((A,R, x0 ), V, T )]](H)(t) =
let G : [0...H) → R

|A|

0

Z

be the solution of G(t ) = x0 +

t0

F (V, T )(G(s))ds
0

(G(t), V, T )
where the above operation reads as follows: ’first line’ = ’third line’, where G is
defined as in ’second line’. If for such an H, G is not unique, then we say that
[[((A, R, x0 ), V, T )]] (H)(t) is ill posed.
In Definition 4 we have explicitly introduced a dependence on a time horizon H,
because it may happen that the solution of the rate equations is defined only for
a finite time horizon [15].
3.1 A Language for Experimental Protocols
Our goal is to build a simple core language that gives an integrated representation
of a discrete protocol within the physical process being implemented on. We
consider the following language modelling the basic operation of a lab protocol.
Definition 5 (Syntax of a Protocol) Given a set of variables V ar, the syntax
of a protocol P for a given fixed CRN C = (A, R) is
P =

x

(sample variable)

(x0 , V, T )

(initial condition)

M ix(P1 , P2 )
let x = P1 in P2
let x, y = Dispense(P1 , p) in P2
Equilibrate(P, t)
Dispose(P )

(mix samples)
(define variable)
(dispense samples)
(let time pass)
(discard P)

where T, V, t ∈ R≥0 , x, y ∈ V ar, p ∈ R(0,1) . Moreover, let-bound variables must
occur exactly once ( that is, be free) in P2 .
A protocol P yields a sample, which is the result of operations of Equilibrate,
Mix, Dispose and Dispense, over a CRS. This syntax allows one to create and
manipulate new samples using Mix (put together different samples), Dispense
(separate samples) and Dispose (discard samples) operations. Note that the CRN
is common for all samples. However, different samples may have different initial
conditions. The single-occurrence (linearity) restriction implies that a sample
cannot be duplicated or eliminated from the pool.

Example 2. We use let x, = Dispense(P1 , p) in P2 as a short-hand for let x, y =
Dispense(P1 , p) in M ix(Dispose(y), x). The protocol (call it P ro1 ) represented
graphically in Figure 1 is defined formally as
P ro1 =let A = ([(H + , 0.1M ); (Cl− , 0.1M )], 1.0mL, 25.0o C) in
let B = ([(N a+ , 0.1M ); (OH − , 0.1M )], 1.0mL, 25.0o C) in
let a, = Dispense(A, p1 ) in
let b, = Dispense(B, p2 ) in
Equilibrate(M ix(a, b), t).
In order to define the semantics of a protocol we introduce the following
definitions.
Definition 6 (Free Variables) The set of Free Variables (FV) of a protocol P
is defined inductively as follows:
F V (x) = {x}
F V ((x0 , V, T )) = {}
F V (M ix(P1 , P2 )) = F V (P1 ) ∪ F V (P2 )
F V (let x = P1 in P2 ) = F V (P1 ) ∪ (F V (P2 ) − {x})
F V (let x, y = Dispense(P1 , p) in P2 ) = F V (P1 ) ∪ (F V (P2 ) − {x, y})
F V (Equilibrate(P, t)) = F V (P )
F V (Dispose(P )) = F V (P ).
We define the operation of substitution of a protocol into a variable as follows.
Definition 7 (Substitution) P1 {x ← P2 } is defined inductively as follows:
x{x ← P } = P
y{x ← P } = y, for x 6= y
M ix(P1 , P2 ){x ← P3 } = M ix(P1 {x ← P3 }, P2 {x ← P3 })
(let x = P1 in P2 ){x ← P3 } = (letx = P1 {x ← P3 } in P2 )
(let y = P1 in P2 ){x ← P3 } = (let y = P1 {x ← P3 }inP2 {x ← P3 }),
for x 6= y and y 6∈ F V (P3 )
Equilibrate(P, t){x ← P1 } = Equilibrate(P {x ← P1 , t})
Dispose(P ){x ← P1 } = Dispose(P {x ← P1 }).
The following equivalences can be shown structurally, based on the definitions
above.
Proposition 1. (Equivalence Relationships)
let x = P1 in P2 = P2 {x ← P1 }
let x = P1 in P2 = let y = P1 in (P2 {x ← y}) for y 6∈ F V (P2 ).

3.2

Deterministic Semantics of a Protocol

We introduce the deterministic semantics for a protocol. Then, in the next Section, we extend such a semantics in order to take into account errors and inaccuracies within the protocol, which in practice may be quite relevant.
The deterministic semantics of a protocol P for a CRN C = (A, R), under
a given environment ρ : V ar → S, is a function [[P ]]ρ : (V ar → S) → S that is
defined inductively as follows.
Definition 8 (Deterministic Semantics of a Protocol) Let S = (R|A| × R≥0 ×
R≥0 ), then the deterministic semantics of a protocol P for CRN C = (A, R),
under environment ρ : V ar → S is defined inductively as follows
[[x]]ρ = ρ(x)
[[x0 , V, T ]]ρ = (x0 , V, T )
[[M ix(P1 , P2 )]]ρ =
let (x10 , V1 , T1 ) = [[P1 ]]ρ
let (x20 , V2 , T2 ) = [[P2 ]]ρ
T1 V1 + T2 V2
x10 V1 + x20 V2
, V1 + V2 ,
)
V1 + V2
V1 + V2
[[let x = P1 in P2 ]]ρ =
(

let (x0 , V, T ) = [[P1 ]]ρ
let ρ1 = ρ{x ← (x0 , V, T )}
[[P2 ]]ρ1
[[let x, y = Dispense(P1 , p) in P2 ]]ρ =
let (x0 , V, T ) = [[P1 ]]ρ
let ρ1 = ρ{x ← (x0 , V · p, T ), y ← (x0 , V · (1 − p), T )}
[[P2 ]]ρ1
[[Equilibrate(P, t)]]ρ =
let (x0 , V, T ) = [[P ]]ρ
[[(A, R, x0 ), V, T )]](H)(t)
[[Dispose(P )]]ρ = (0|Λ| , 0, 0),
where H ∈ R≥0 is such that for any Equilibrate(P, t), [[(A, R), x0 , V, T )]](H)(t)
is well posed. If such an H does not exist, we say that P is ill posed.
The above semantics identifies a protocol which outputs the concentration of the
species, the volume, and the temperature of the sample at final time.
3.3

Deterministic Semantics of a Protocol as a PDMP

Given a protocol P and an environment ρ, [[P ]]ρ induces semantics that correspond to the solution of a PDMP H = (Q, d, G, F, Λ, R) as per Definitions

1 and 2. In the corresponding PDMP model H, Q represents the set of discrete operations, so that d(q) = |A| + 1 denotes the continuous dimension (the
number of continuous variables). The vector field F is given by Definition 4,
with an additional clock variable time representing time as dtime
= 1. For each
dt
Equilibrate(P, t) step, there is a guard set defined as time ≥ t: when a trajectory
enters this guard, an associated reset is modeled with the identity function. The
resulting process is a PDMP without random jumps (that is, where Λ(q, x) = 0
for any q ∈ Q, x ∈ R|A|+1 ) and with non-probabilistic resets R. As such, the
resulting PDMP is also a (non-probabilistic) hybrid model [2]. We elucidate this
with the next example.
Example 3. Consider the protocol P ro1 introduced in Example 2. The CRN of
the system comprises the reactions given in the CRN in Example 1. According
to Definition 11, the state of variable H + in time is given by the solution of the
following equation:
Z t
+
+
H (t) = H (0) −
kN a+ (s)OH − (s)H + (s)Cl− (s)ds,
0
−7.4

2 10
. Then, given an environment ρ, [[P ro1 ]]ρ is the
where H + (0) = p1 0.1+p
p1 +p2
solution of the PDMP H = (Q, d, G, F, Λ, R), where Q = {q}; d(q) = R8 ;
G(q, x) = 1 iff xtime ≥ t, where xtime is the component of x relative to variable time; and F is given by Definition 11. Notice again that Λ(q, x) = 0 for any
x and R(q, x) = δ(x), where δ(x) is the Dirac delta function centered at x.

3.4

Stochastic Semantics of a Protocol, and Interpretation as a
PDMP
The semantics of Definition 11 are fully deterministic, and indeed are shown to
map into a fully non-probabilistic PDMP model. However, it is often the case
that operations of Dispense and Equilibrate are stochastic in nature, due to
the fact that they are performed by humans and in view of experimental inaccuracies related to lab equipment. In what follows, we encompass this features by
extending the previously defined semantics with stochasticity. More precisely,
– in the Equilibrate(P, t) step, time is sampled from a distribution;
– the resulting volume after a Dispense step is sampled from a distribution.
The first characteristic models the fact that in real experiments the system is not
equilibrated for exactly t seconds, as it may start or be stopped at different time
instants, and accounts for the fact that after a mix of samples well mixed conditions are not reached instantaneously, whereas the second takes into account the
error of pipetting devices, whose ranges and parameters have been standardized
t0
(standard ISO 8655). For the first feature, consider the function T (t0 , t) = e− t ,
defined for two values t0 , t ∈ R≥0 . This function corresponds to the density
function of an exponential random variable, modelling random arrivals. For the
m
second function, let B(Rm
≥0 ) be the Borel sigma-algebra over R≥0 , m > 0. Then,

we consider the following function D : B(R(0,1) ) × R≥0 × R[0,1] → [0, 1], which
assigns to D(·, V, p) a probability measure in B(R[0,1] ). Function D is used to
reset the volume randomly, after a discrete operation. (As an anticipation of
results, notice that both functions T and D can be mapped to elements in the
syntax in a PDMP model.)
We define the Stochastic Semantics of a protocol as an extension of the
deterministic ones in Definition 11. For the sake of compactness, we explicitly
write only the operators that differ from the earlier ones.
Definition 9 (Stochastic Semantics of a Protocol) Let S = (R|A| × R≥0 × R≥0 ),
then the semantics of a protocol P for CRN C = (A, R), under environment
ρ : V ar → S and functions T , D, as defined above, is defined inductively as
follows
[[let x, y = Dispense(P1 , p) in P2 ]]ρ =
let (x0 , V, T ) = [[P1 ]]ρ
let p0 being sampled f rom D(·, V, p)
let ρ1 = ρ{x ← (x0 , V · p0 , T ), y ← (x0 , V · (1 − p0 ), T )}
[[P2 ]]ρ1
[[Equilibrate(P, t)]]ρ =
let (x0 , V, T ) = [[P ]]ρ
let I be a R≥0 − valued random variable such that for s ∈ R≥0
P rob(I > s) = T (s, t)
[[(A, R, x0 ), V, T )]](H)(I),
where H ∈ R≥0 is such that for any Equilibrate(P, t), and any I random variable
such that P rob(I > s) = T (s, t), [[(A, R), x0 , V, T )]](H)(I) is well posed with
probability 1. If such an H does not exist, we say that [[P ]]ρ is ill posed.
D is a transition kernel that depends only on the current state of the system. T
is the cumulative probability distribution of a random variable with exponential
distribution. As a consequence, according to Definition 2, [[P ]]ρ induces semantics
that are again solution of a PDMP. However, here T determines the probability
of changing discrete state and D acts as a probabilistic reset, there are no guards,
and the continuous dynamics evolve according to the ODE in Definition 4.
Next, we leverage results from the analysis of PDMP models and export
them over the protocol language. The following assumptions guarantee that the
solution of the PDMP induced by [[P ]]ρ exists, establish that it is a strong Markov
process, and allow to exclude pathological Zeno behaviours [12,19].
Assumption 1
– Let A0 , A1 ⊂ B(R[0,1] ) be the smallest sets in B(R[0,1] ) containing respectively
0 and 1. Then, D(A0 , V, p) = D(A1 , V, p) = 0 for any p ∈ (0, 1), V 6= 0. That
is, the Volume after a dispense is zero with probability zero.

– Let F : R|A| → R|A| be the drift term of the rate equations (Eqn (4)). Then,
F is a globally Lipschitz function.
– For any Equilibrate(·, t) we have that t > 0.
Let us interpret these assumptions over the protocol languages. The first assumption guarantees that the volume of a non-empty sample cannot be 0 almostsurely. The second assumption guarantees that the solution of (4) exists and
does not hit infinity in finite time. This excludes non-physical reactions like
X + X → X + X + X. The third assumption guarantees that we have a finite
number of jumps over a finite time, thus excluding Zeno behaviours [12,13].
Example 4. Consider the protocol introduced in Example 2. For σ1 > 0, A ⊂
−

R[0.1,0.8] . Assume that D(A, p, V̄ ) =

R
A

e

R 0.8
0.1

x−p
2
2σ1

−

e

dx

x−p
2
2σ1

. That is, D(·, p, V ) is a trun-

dx

cated Gaussian measure centered at p and independent of the volume. Then,
according to Definition 9, we have the following Stochastic Semantics:
Z I
H + (I) =H + (0) −
kN a+ (s)OH − (s)H + (s)Cl− (s)ds,
0
−7.4

2 10
. Here I is a random variable with an exponential
with HCl(0) = V1 0.1+V
V1 +V2
distribution with rate T1 , V1 is a random variable sampled from D(·, p1 , 1), and
V2 is a random variable sampled from D(·, p2 , 1).

4

Extending the Protocol Language with Observations

The language introduced in Section 3 can be extended in a number of directions,
according to specific envisioned scenarios for the protocols. A common task is to
take observations of the state of the protocol. That is, often it is useful to store
the state of the system at different times or when a particular event happens.
As some of the events may be stochastic, in general, it is not possible to know
before the simulation starts when a particular event happens. Consequently,
observations need to be included in the language.
Definition 10 (Extended Syntax). Given a set of variables V ar, the syntax of
a protocol for a given fixed CRN C = (A, R) and idn ∈ N is
P =

x

(sample variable)

(x0 , V, T )

(initial condition)

M ix(P1 , P2 )
let x = P1 in P2
let x, y = Dispense(P1 , p) in P2
Equilibrate(P, t)
Dispose(P )
Observe(P, idn)

(mix samples)
(define variable)
(dispense samples)
(let time pass)
(discard P)
(observe sample)

where T, V, t ∈ R≥0 , x, y ∈ V ar, p ∈ R[0,1] . Moreover, let-bound variables
must occur exactly once ( that is, be free) in P2 .

Observe(P, idn) makes an observation of protocol P after its execution, and
identifies such an observation with identifier idn. In order to include observations
we extend the semantics as detailed next, where we consider in detail just the
deterministic semantics, focusing on a few key operators. The other operators
and the Stochastic Semantics follow similarly.
Definition 11 (Extended Deterministic Semantics) For CRN C = (A, R) let
S = R|A| × R≥0 × R≥0 , Obs = R≥0 × N × R|A| , Obs∗ , an eventually empty
set of Obs and M = S × Obs∗ × R≥0 . The semantics of a protocol P , under
environment ρ : V ar → M, is a function [[P ]] : (V ar → M) × R≥0 → M defined
inductively as follows
[[M ix(P1 , P2 )]]ρt =
let ((x10 , V1 , T1 ), Obs1 , t1 ) = [[P1 ]]ρt
let ((x20 , V2 , T2 ), Obs2 , t2 ) = [[P2 ]]ρt
T1 V1 + T2 V2
x10 V1 + x20 V2
, V1 + V2 ,
), Obs1 :: Obs2 , max(t1 , t2 ))
V1 + V2
V1 + V2
[[Observe(P, idn)]]ρt =
((

let ((x0 , V, T ), Obs, t1 ) = [[P ]]ρt
let O = (x0 , idn, t1 )
((x0 , V, T ), Obs ∪ O, t1 )
[[Equilibrate(P, t)]]ρt0 =
let ((x0 , V, T ), Obs, t1 ) = [[P ]]ρt0
([[(A, R), x0 , V, T )]](H)(t), Obs, t1 + t),
where H ∈ R≥0 is such that for any Equilibrate(P, t), [[(A, R), x0 , V, T )]](H)(t)
is well posed. If such H does not exist, we say that P is ill posed.
Note that the above syntax does not prevent the programmer to assign the same
identifier to two distinct observations. We further stress that often observations
of the state of an experiment are not exact, but corrupted by sensing noise.
For instance, this is what happens with noisy fluorescence measurements. This
noise can be easily taken into account at a semantical level by sampling an
observation from a distribution with added noise, where the noise level depends
on the particular measure technique or instrumentation. Finally, we can as well
extend the sample semantic to take into account noise in Dispense operations.

5

Case Study

As a case study we consider an experimental protocol for DNA strand displacement. DNA strand displacement (DSD) is a design paradigm for DNA
nano-devices [10]. In such a paradigm, single-stranded DNA acts as signals and
double-stranded (or more complex) DNA structures act as gates. The interactions between signals and gates allow one to generate computational mechanisms
that can operate autonomously at the molecular level [26]. The DSD programming language has been developed as a means of formally programming and

analyzing such devices [20,10]. Here, we consider an AN D circuit implemented
in DSD, which can be represented with the reactions in Figure 2b. Strands
Input1 =< 1∗ 2 > and Input2 =< 3 4∗ > represent the two inputs, while
strand Output =< 2 3 > is the output. Strand Gate = {1∗ }[2 3]{4∗ } is an
auxiliary strand. The Output strand is released only if both the inputs react
with the Gate gate. We consider the protocol in Figure 2a, which can be written
formally as follows. We use let x,= Dispense(P 1, p) in P 2 as a short-hand for
let x, y = Dispense(P 1, p) in M ix(Dispose(y), x)
P1 =let In1 = ((Input1, 100.0nM ), 0.1mL, 25.0o C) in
let In2 = ((Input2, 100.0nM ), 0.1mL, 25.0o C) in
let GA = ((Output, 100.0nM ), 0.1mL, 25.0o C) in
let GB = ((GateB , 100.0nM ), 0.1mL, 25.0o C) in
let sGA,= Dispense(GA, p1 ) in
let sGB,= Dispense(GB, p2 ) in
let sIn1,= Dispense(In1, p3 ) in
let sIn2,= Dispense(In1, p4 ) in
Observe(Equilibrate(M ix(M ix(Equilibrate(M ix(sGA, sGB), t1 ), sIn1), sIn2), t2 ), idn).
The protocol proceeds as follow: Output and GateB strands are dispensed from
the original samples. Then, they are let evolve for t1 seconds to create Gate
strands. Then, the two inputs are dispensed from their samples. The resulting
samples are mixed and the resulting solution evolves for t2 seconds. Finally, we
collect the final sample, observe the results, and associate to the observation
the identifier ’idn’. According to the standard ISO 8655 for a volume of 1mL,

Fig. 2: (A) Graphical representation of the protocol. (B) Graphical representation
of the reactions between the different DNA strands in the considered solution.
For example, in the second reaction, stand {1∗ }[2 3]{4∗ } reacts with < 1∗ 2 >
at a rate 0.0003, and there exists an inverse reaction with rate 0.1126.
the maximum standard deviation of a particular pipetting device is 0.3µL per

single operation. To incorporate such an error in our model, we make use of the
Stochastic Semantics. Thus, the concentration of the Output strand at the end
of the protocol is a random variable. It is common that the reaction rates of the
physical system are not known exactly and they may be affected by extrinsic
noise [23]. This leads to another source of uncertanity in the output of the
protocol. To estimate the distribution of the output under the effect of both
these sources of noise we extend our semantics to sample the rate of a reaction
from a normal distribution with variance equals to half of its mean (sub-Poisson
noise). In Figure 3a we plot 4500 executions resulting from the protocol. From
the figure it is easy to realize how the difference source of noise may have a
distinctive effect on the final outcome of the experiments.
In many experimental protocols, one of the key challenges is to synthesize
their optimal discrete parameters, to optimize the probability of obtaining desired behaviours. Here, we assume perfect knowledge of the reaction rates of the
physical system and p1 = p2 = 0.4. Our goal is to see how the concentration of
the Output changes varying (p3 , p4 ) ∈ [0.45, 0.65]×[0.45, 0.65]. We are interested
in the following property
PSaf e ([3.0·10−4 , 3.5·10−4 ]) = P rob(Output(t0 ) ∈ [3.0·10−4 , 3.5·10−4 ]|t0 = tf inal ),
where tf inal is the final time of the protocol. The following probability is estimated using Statistical Model Checking [22] in Figure 3b, which in this context
reduces to Monte-Carlo sampling. From Figure 3b it is easy to infer that the

Fig. 3: (A): (red) 1500 execution of the protocol assuming the physical model
is fully known, and the only source of noise is the discrete parameters of the
protocols (p1 , p2 , p3 , p4 ). (yellow) 1500 executions of the protocol when the rates
of the physical system are sampled from a sub-Poisson distribution, and discrete
operations are exact. (blue) 1500 simulations of the protocol when both sources
of noise are active.(B): PSaf e ([3.0 · 10−4 , 3.5 · 10−4 ]) as a function of p3 and p4 .
Each cell is estimated from 20000 executions of the protocol.
optimal value for such property is not unique (it is attained at values over the
yellow band) and obtained, for instance, at (p3 , p4 ) = (0.5, 0.54).

6

Discussion

We presented a language to model experimental biological protocols, and provided semantics to this protocol in terms of PDMPs. That is, to each experimental biological protocol is associated a particular instance of a stochastic hybrid

process. Our language provides a unified description of the model of the system
being experimented on, together with the discrete events representing the parts
of biological protocols dealing with handling samples. Moreover, we allow the
modeller to take into account uncertainties in both the model structure and the
equipment tolerances. This makes our language a suitable tool for both experimental and computational biologists. Our objective has been that of providing
a basic language, yielding an integrated representation of an experimental biological protocol. To this end, we have kept the language as simple as possible,
showing how different extensions are easy to be integrated. For instance, in our
denotational semantics the dynamics of a physical process is given by a set of
ODEs. This is accurate when the number of molecules involved is big enough,
as in the discussed example of DNA strand displacement (DSD). However, in
other scenarios, such as localized computation or gene expression, this might be
unsatisfactory as stochasticity becomes important [5,14]: nevertheless, the semantics presented here can be easily extended to incorporate such stochasticity,
which can be done for example by considering more general classes of stochastic
hybrid processes, such as switching diffusions [21,27]. Another relatively simple
extension is to include finite loops or operations based on concentrations.
One of the main advantages in providing a language with formal semantics for
experimental protocols is that experimental protocols can now be quantitatively
analyzed inexpensively in-silico, and classical problems of analysis of CRNs, such
as parameter estimation [6], can be performed within this framework, also by
taking into account the discrete operations of the protocol, which influences the
dynamics of the system. An additional target is to provide automated techniques
to synthesise optimal protocols, or protocols that are guaranteed to perform as
desired. This can be attained by tapping into the mature literature on formal
verification and strategy synthesis of PDMPs, or that of any other more special model that the given protocol can be mapped onto. Notions of finite-state
abstractions [28] and of probabilistic bisimulations [1], as well as algorithms for
probabilistic model checking of stochastic hybrid models [21] will be relevant
towards this goal.

References
1. A. Abate. Probabilistic bisimulations of switching and resetting diffusions. In
Proceedings of the 49th IEEE Conference of Decision and Control, pages 5918–
5923, 2010.
2. R. Alur, C. Courcoubetis, N. Halbwachs, T. A. Henzinger, P.-H. Ho, X. Nicollin,
A. Olivero, J. Sifakis, and S. Yovine. The algorithmic analysis of hybrid systems.
Theoretical computer science, 138(1):3–34, 1995.
3. V. Ananthanarayanan and W. Thies. Biocoder: A programming language for standardizing and automating biology protocols. Journal of Biological Engineering,
4(1):13, Nov 2010.
4. A. Andreychenko, L. Mikeev, D. Spieler, and V. Wolf. Parameter identification for
markov models of biochemical reactions. In Computer Aided Verification, pages
83–98. Springer, 2011.

5. A. Arkin, J. Ross, and H. H. McAdams. Stochastic kinetic analysis of developmental pathway bifurcation in phage λ-infected Escherichia coli cells. Genetics,
149(4):1633–1648, 1998.
6. L. Cardelli, M. Češka, M. Fränzle, M. Kwiatkowska, L. Laurenti, N. Paoletti, and
M. Whitby. Syntax-guided optimal synthesis for chemical reaction networks. In International Conference on Computer Aided Verification, pages 375–395. Springer,
2017.
7. L. Cardelli, M. Kwiatkowska, and L. Laurenti. Programming discrete distributions
with chemical reaction networks. In International Conference on DNA-Based Computers, pages 35–51. Springer, 2016.
8. L. Cardelli, M. Kwiatkowska, and L. Laurenti. Stochastic analysis of chemical
reaction networks using linear noise approximation. Biosystems, 149:26–33, 2016.
9. L. Cardelli, M. Kwiatkowska, and L. Laurenti. A stochastic hybrid approximation
for chemical kinetics based on the linear noise approximation. In International Conference on Computational Methods in Systems Biology, pages 147–167. Springer,
2016.
10. Y.-J. Chen, N. Dalchau, N. Srinivas, A. Phillips, L. Cardelli, D. Soloveichik, and
G. Seelig. Programmable chemical controllers made from DNA. Nature nanotechnology, 8(10):755–762, 2013.
11. N. Dalchau, N. Murphy, R. Petersen, and B. Yordanov. Synthesizing and tuning
chemical reaction networks with specified behaviours. In International Workshop
on DNA-Based Computers, pages 16–33. Springer, 2015.
12. M. H. Davis. Piecewise-deterministic Markov processes: A general class of nondiffusion stochastic models. Journal of the Royal Statistical Society. Series B
(Methodological), pages 353–388, 1984.
13. M. H. Davis. Markov Models & Optimization, volume 49. CRC Press, 1993.
14. K. E. Dunn, F. Dannenberg, T. E. Ouldridge, M. Kwiatkowska, A. J. Turberfield,
and J. Bath. Guiding the folding pathway of dna origami. Nature, 525(7567):82,
2015.
15. S. N. Ethier and T. G. Kurtz. Markov processes: characterization and convergence,
volume 282. John Wiley & Sons, 2009.
16. M. Feinberg. Chemical reaction network structure and the stability of complex
isothermal reactors – The deficiency zero and deficiency one theorems. Chemical
Engineering Science, 42(10):2229–2268, 1987.
17. L. P. Freedman, I. M. Cockburn, and T. S. Simcoe. The economics of reproducibility in preclinical research. PLOS Biology, 13(6):1–9, 06 2015.
18. E. C. Hayden. The automated lab. Nature, 516(7529):131–132, 12 2014.
19. P. Kouretas, K. Koutroumpas, J. Lygeros, and Z. Lygerou. Stochastic hybrid
modeling of biochemical processes. Stochastic Hybrid Systems, 24(9083), 2006.
20. M. R. Lakin, S. Youssef, F. Polo, S. Emmott, and A. Phillips. Visual DSD: a
design and analysis tool for DNA strand displacement systems. Bioinformatics,
27(22):3211–3213, 2011.
21. L. Laurenti, A. Abate, L. Bortolussi, L. Cardelli, M. Ceska, and M. Kwiatkowska.
Reachability computation for switching diffusions: Finite abstractions with certifiable and tuneable precision. In Proceedings of the 20th International Conference
on Hybrid Systems: Computation and Control, pages 55–64. ACM, 2017.
22. A. Legay, B. Delahaye, and S. Bensalem. Statistical model checking: An overview.
RV, 10:122–135, 2010.
23. J. Paulsson. Summing up the noise in gene networks. Nature, 427(6973):415–418,
2004.

24. M. I. Sadowski, C. Grant, and T. S. Fell. Harnessing qbd, programming languages,
and automation for reproducible biology. Trends in Biotechnology, 34(3):214–227,
10 2017.
25. D. S. Scott and C. Strachey. Toward a mathematical semantics for computer
languages, volume 1. Oxford University Computing Laboratory, Programming
Research Group, 1971.
26. G. Seelig, D. Soloveichik, D. Y. Zhang, and E. Winfree. Enzyme-free nucleic acid
logic circuits. science, 314(5805):1585–1588, 2006.
27. G. Yin and C. Zhu. Hybrid switching diffusions: properties and applications, volume 63. Springer New York, 2010.
28. M. Zamani and A. Abate. Symbolic models for randomly switched stochastic
systems. Systems and Control Letters, 69:38–46, 2014.

