A Language for Modeling and Optimizing Experimental Biological Protocols

Abstract

Automation is becoming ubiquitous in all laboratory activities, moving towards precisely defined and codified laboratory protocols. However, the integration between laboratory protocols and mathematical models is still lacking. Models describe physical processes, while protocols define the steps carried out during an experiment: neither cover the domain of the other, although they both attempt to characterize the same phenomena. We should ideally start from an integrated description of both the model and the steps carried out to test it, to concurrently analyze uncertainties in model parameters, equipment tolerances, and data collection. To this end, we present a language to model and optimize experimental biochemical protocols that facilitates such an integrated description, and that can be combined with experimental data. We provide probabilistic semantics for our language in terms of Gaussian processes (GPs) based on the linear noise approximation (LNA) that formally characterizes the uncertainties in the data collection, the underlying model, and the protocol operations. In a set of case studies, we illustrate how the resulting framework allows for automated analysis and optimization of experimental protocols, including Gibson assembly protocols.

1. Introduction

Automation is becoming ubiquitous in all laboratory activities: protocols are run under reproducible and auditable software control, data are collected by high-throughput machinery, experiments are automatically analyzed, and further experiments are selected to maximize knowledge acquisition. However, while progress is being made towards the routine integration of sophisticated end-to-end laboratory workflows and towards the remote access to laboratory facilities and procedures [1,2,3,4,5], the integration between laboratory protocols and mathematical models is still lacking. Models describe physical processes, either mechanistically or by inference from data, while protocols define the steps carried out during an experiment in order to obtain experimental data. Neither models nor protocols cover the domain of the other, although they both attempt to characterize the same phenomena. As a consequence, it is often hard to attribute causes of experimental failures: whether an experiment failed because of a misconception in the model, or because of a misstep in the protocol. To confront this problem, we need an approach that integrates and accounts for all the components, theoretical and practical, of a laboratory experiment. We should ideally start from an integrated description from which we can extract 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. This is essential to enable automated model synthesis and falsification by concurrently taking into account uncertainties in model parameters, equipment tolerances, and data collection.

We present a language to model and optimize experimental biochemical protocols that provides such an integrated description of the protocol and of the underlying molecular process, and that can be combined with experimental data. From this integrated representation, both the model of a phenomenon and the steps carried out to test it can be separately extracted. Our approach is illustrated in Figure 1.

We provide probabilistic semantics for our language in terms of a Gaussian process (GP) [6], which can be used to characterize uncertainties. Such a semantics arises from a linear noise approximation (LNA) [7,8] of the dynamics of the underlying biochemical system and of the protocol operations, and it corresponds to a Gaussian noise assumption. We show that in a Bayesian framework, the resulting semantics can be combined in a principled way with experimental data to build a posterior process that integrates our prior knowledge with new data. We demonstrate that the Gaussian nature of the resulting process allows one to efficiently and automatically optimize the parameters of an experimental protocol in order to maximize the performances of the experiment. On a series of case studies, including a Gibson Assembly protocol [9] and a Split and MIx protocol, we highlight the usefulness of our approach and how it can have an impact on scientific discovery.

Related work and novelty.

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 [10]. Second, the emergence of ‘cloud labs’ 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 BioScript [11,12], BioCoder [2], Autoprotocol [13], and Antha [14]. These languages provide expressive, high-level protocol descriptions but consider each experimental sample as a labeled ‘black-box’. This makes challenging the 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 we capture the details of experimental samples, enabling us to track properties of chemical solutions (concentrations and temperatures) as the chemicals that react during the execution of a protocol. This allows us to formalize and verify requirements for the correct execution of a protocol and to optimize various protocol or system parameters to satisfy these specifications.

Our language derives from our previous conference paper [5] by introducing stochasticity in chemical evolution, and by providing a Gaussian semantics, particularly for protocol operations, for the effective analysis of the combined evolution of stochastic chemical kinetics and protocols and experimental data. The language semantics has been already implemented (but not previously formally presented) in a chemical/protocol simulator [15]. We expect that this style of semantics could be adapted to other protocol languages, along the principles we illustrate, to provide a foundation for the analysis of complex stochastic protocols in other settings.

Paper outline.

In what follows, we first introduce the syntax and semantics of our language. We then show how we can perform optimization of the resulting Gaussian process integrated with data. Finally, we illustrate the usefulness of our approach on several case studies, highlighting the potential impact of our work on scientific discovery.

2. Materials and Methods

In this section, we introduce the syntax of the language we propose for modeling experimental protocols. A formal semantics of the language, based on denotational semantics [16], is then discussed. The physical process underlying a biological experimental protocol is modeled as a chemical reaction system (CRS). As a consequence, before introducing the language for experimental protocols, we first formally introduce chemical reaction networks (CRNs) and chemical reaction systems (CRSs).

2.1. Chemical Reaction Network (CRN)

Chemical reaction networks (CRNs) is a standard language for modeling biomolecular interactions [17].

Definition 1

(CRN). A chemical reaction network $\mathcal{C}=(\mathcal{A},\mathcal{R})\in CRN=\Lambda \times \mathcal{T}$ is a pair of a finite set $\mathcal{A}\subseteq \Lambda$ of chemical species, of size $\left|\mathcal{A}\right|$, and a finite set $\mathcal{R}\subseteq \mathcal{T}$ of chemical reactions. A reaction $\tau \in \mathcal{R}$ is a triple $\tau =(r_{\tau },p_{\tau },k_{\tau })$, where $r_{\tau }\in {\mathbb{N}}^{\left|\Lambda \right|}$ is the source complex, $p_{\tau }\in {\mathbb{N}}^{\left|\Lambda \right|}$ is the product complex and $k_{\tau }\in {\mathbb{R}}_{>0}$ is the coefficient associated with the rate of the reaction. The quantities $r_{\tau }$ and $p_{\tau }$ are the stoichiometry of reactants and products. The stoichiometric vector associated to τ is defined by ${\upsilon }_{\tau }=p_{\tau }-r_{\tau }$. Given a reaction ${\tau }_i=([1,0,1],[0,2,0],k_i)$ we refer to it visually as ${\tau }_i:{\lambda }_1+{\lambda }_3{\to }^{k_i}2{\lambda }_2$.

Definition 2

(CRN State). Let $\mathcal{C}=(\mathcal{A},\mathcal{R})\in CRN$. A state of $\mathcal{C}$, of the form $(\mu ,\Sigma ,V,T)$$\in S=({\mathbb{R}}_{\ge 0}^{\left|\mathcal{A}\right|}\times {\mathbb{R}}^{\left|\mathcal{A}\right|\times \left|\mathcal{A}\right|}\times {\mathbb{R}}_{\ge 0}\times {\mathbb{R}}_{\ge 0})$, consists of a concentration vector μ (moles per liter), a covariance matrix Σ, a volume (liters), and a temperature (degrees Celsius). A CRN together with a (possibly initial) state is called a Chemical Reaction System (CRS) having the form $(\mathcal{A},\mathcal{R}),(\mu ,\Sigma ,V,T)$.

A Gaussian process is a collection of random variables, such that every finite linear combination of them is normally distributed. Given a state of a CRN, its time evolution from that state can be described by a Gaussian process indexed by time whose mean and variance are given by a linear noise approximation (LNA) of the Chemical Master Equation (CME) [7,18] and is formally defined in the following definitions of CRN Flux and CRN Time Evolution.

Definition 3

(CRN Flux). Let $(\mathcal{A},\mathcal{R})$ be a CRN. Let $F(V,T)\in {\mathbb{R}}_{\ge 0}^{\left|\mathcal{A}\right|}\to {\mathbb{R}}^{\left|\mathcal{A}\right|}$ be the flux of the CRN at volume $V\in {\mathbb{R}}_{\ge 0}$ and temperature $T\in {\mathbb{R}}_{\ge 0}$. For a concentration vector $\mu \in {\mathbb{R}}_{\ge 0}^{\left|\mathcal{A}\right|}$ we assume $F(V,T)\left(\mu \right)={\sum }_{\tau \in \mathcal{R}}{\upsilon }_{\tau }{\alpha }_{\tau }(V,T,\mu )$, with stoichiometric vector ${\upsilon }_{\tau }$ and rate function ${\alpha }_{\tau }$. We call $J_F$ the Jacobian of $F(V,T)$, and $J_F^{\top }$ its transpose. Further, define $W(V,T)\left(\mu \right)={\sum }_{\tau \in \mathcal{R}}{\upsilon }_{\tau }{\upsilon }_{\tau }^{\top }{\alpha }_{\tau }(V,T,\mu )$ to be the diffusion term.

We should remark that, if one assumes mass action kinetics, then for a reaction $\tau$ it holds that ${\alpha }_{\tau }(V,T,\mu )=k_{\tau }^{\prime }(k_{\tau },V,T){\mu }^{r_{\tau }}$ where ${\mu }^s={\prod }_{i=1}^{\left|\Lambda \right|}{\mu }_i^{s_i}$ is the product of the reagent concentrations, and $k_{\tau }^{\prime }$ encodes any dependency of the reaction rate $k_{\tau }$ on volume and temperature. Then, the Jacobian evaluated at $\mu$ is given by $J_{F_{ik}}\left(\mu \right)=\frac{\partial F(V,T){\left(\mathbf{x}\right)}_i}{\partial {\mathbf{x}}_k}{|}_{\mathbf{x}=\mu }={\sum }_{\tau \in \mathcal{R}}{\upsilon }_{{\tau }_i}{r}_{{\tau }_k}\frac{{\alpha }_{\tau }(V,T,\mu )}{{\mu }_k}$ for $i,k\in \left|\mathcal{A}\right|$ [8].

In the following, we use $\mu ,\Sigma$ for concentration mean vectors and covariance matrices, respectively, and $\mathbf{\mu },\Sigma$ (boldface) for the corresponding functions of time providing the evolution of means and covariances.

Definition 4

(CRN Time Evolution). Given a CRS $(\mathcal{A},\mathcal{R}),(\mu ,\Sigma ,V,T)$, its evolution at time $t<H$ (where $H\in {\mathbb{R}}_{\ge 0}\cup \{\infty \}$ is a time horizon) is the state $({\mathbf{\mu }}_{\mu }\left(t\right),{\Sigma }_{\mu ,\Sigma }\left(t\right),V,T)$ obtained by integrating its flux up to time t, where:

$$\begin{array}{cc}\hfill & {\mathbf{\mu }}_{\mu }\left(t\right)=\mu +{\int }_0^{t}F(V,T)\left({\mathbf{\mu }}_{\mu }\left(s\right)\right)ds\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & {\Sigma }_{\mu ,\Sigma }\left(t\right)=\Sigma +{\int }_0^t{J}_F\left({\mathbf{\mu }}_{\mu }\left(s\right)\right){\Sigma }_{\mu ,\Sigma }\left(s\right)+{\Sigma }_{\mu ,\Sigma }\left(s\right)J_F^{\top }\left({\mathbf{\mu }}_{\mu }\left(s\right)\right)+W(V,T)\left({\mathbf{\mu }}_{\mu }\left(s\right)\right)ds,\hfill \end{array}$$

(2)

with ${\mathbf{\mu }}_{\mu }\left(0\right)=\mu$ and ${\Sigma }_{\mu ,\Sigma }\left(0\right)=\Sigma$. If, for such an H, $\mathbf{\mu }$ or Σ are not unique, then we say that the evolution is ill-posed. Otherwise, ${\mathbf{\mu }}_{\mu }\left(t\right)$ and ${\Sigma }_{\mu ,\Sigma }\left(t\right)$ define a Gaussian process with that mean and covariance matrix for $t<H$.

An ill-posed problem may result from technically expressible but anomalous CRS kinetics that does not reflect physical phenomena, such as deterministic trajectories that are not uniquely determined by initial conditions, or trajectories that reach infinite concentrations in finite time.

Example 1.

Consider the CRN $\mathcal{C}=\left(\right\{a,b,c\},\mathcal{R})$, where $\mathcal{R}=\left\{\right([1,1,0],[0,2,0],0.2),([0,1,1],$$[0,0,2],0.2\left)\right\}$. Equivalently, $\mathcal{R}$ can be expressed as

$$a+b{\to }^{0.2}b+bb+c{\to }^{0.2}c+c.$$

Then, we have

$$\begin{array}{cc}\hfill & F(a,b,c)=0.2\left[\begin{array}{c}-a·b\\ a·b-b·c\\ b·c\end{array}\right]\hfill \\ \hfill & W(a,b,c)=0.2\left[\begin{array}{ccc}a·b& -a·b& 0\\ -a·b& a·b+b·c& -b·c\\ 0& -b·c& b·c\end{array}\right]\hfill \\ \hfill & J_F(a,b,c)=0.2\left[\begin{array}{ccc}-b& -a& 0\\ b& a-c& -b\\ 0& c& b\end{array}\right]\hfill \end{array}$$

Consider an initial condition of ${\mu }_a=0.1,$${\mu }_b={\mu }_c=0.001$, then the time evolution of mean and variance for the GP defined in Definition 4 are reported in Figure 2.

2.2. A Language for Experimental Biological Protocols

In what follows, we introduce the syntax of our language for experimental biological protocols.

Definition 5

(Syntax of a Protocol). Given a set of sample variables $\mathit{x}\in \mathit{V}\mathit{a}\mathit{r}$ and a set of parameter variables $\mathit{z}\in \mathit{P}\mathit{a}\mathit{r}$, the syntax of a protocol $\mathit{P}\in \mathit{P}\mathit{r}\mathit{o}\mathit{t}$ for a given fixed CRN $\mathcal{C}=(\mathcal{A},\mathcal{R})$ is -4.6cm0cm

$$\begin{array}{cc}\hfill \mathit{P}=\mathit{x}& (\mathit{a}\mathit{sample}\mathit{variable})\hfill \\ \hfill ({\mathit{p}}_1\dots {\mathit{p}}_{\left|\mathcal{A}\right|},{\mathit{r}}_{\mathit{V}},{\mathit{r}}_{\mathit{T}})& (\mathit{a}\mathit{sample}\mathit{with}\mathit{initial}\mathit{concentrations},\mathit{volume},\mathit{temperature})\hfill \\ \hfill \mathit{l}\mathit{e}\mathit{t}\mathit{x}={\mathit{P}}_1\mathit{i}\mathit{n}{\mathit{P}}_2& (\mathit{introduce}\mathit{a}\mathit{local}\mathit{sample}\mathit{variable}\mathit{x})\hfill \\ \hfill \mathit{M}\mathit{i}\mathit{x}({\mathit{P}}_1,{\mathit{P}}_2)& (\mathit{mix}\mathit{samples})\hfill \\ \hfill \mathit{l}\mathit{e}\mathit{t}{\mathit{x}}_1,{\mathit{x}}_2=\mathit{S}\mathit{p}\mathit{l}\mathit{i}\mathit{t}({\mathit{P}}_1,\mathit{p})\mathit{i}\mathit{n}{\mathit{P}}_2& (\mathit{split}\mathit{a}\mathit{sample}{\mathit{P}}_1\mathit{by}\mathit{a}\mathit{proportion}\mathit{p}\mathit{in}(0\cdots 1))\hfill \\ \hfill \mathit{E}\mathit{q}\mathit{u}\mathit{i}\mathit{l}\mathit{i}\mathit{b}\mathit{r}\mathit{a}\mathit{t}\mathit{e}(\mathit{P},\mathit{p})& (\mathit{equilibrate}\mathit{a}\mathit{sample}\mathit{for}\mathit{p}\mathit{sec}\mathit{onds})\hfill \\ \hfill \mathit{D}\mathit{i}\mathit{s}\mathit{p}\mathit{o}\mathit{s}\mathit{e}\left(\mathit{P}\right)& (\mathit{discard}\mathit{sample})\hfill \\ \\ \hfill \mathit{p}=\mathit{z}& (\mathit{a}\mathit{parameter}\mathit{variable})\hfill \\ \hfill \mathit{r}& (\mathit{a}\mathit{literal}\mathit{non}\text{-}\mathit{negative}\mathit{real}\mathit{number})\hfill \end{array}$$

Moreover, let-bound variables x, $x_1$, $x_2$ must occur (as free variables) exactly once in $P_2$.

A protocol P manipulates samples (which are CRN states as in Definition 2) through a set of operations, and finally yields a sample as the result. This syntax allows one to create and manipulate new samples using Mix (put together different samples), Split (separate samples) and Dispose (discard samples) operations. Note that the CRN is common to all samples, but different samples may have different initial conditions and hence different active reactions. The single-occurrence (linearity) restriction of sample variables implies that a sample cannot be duplicated or forgotten.

2.3. Gaussian Semantics for Protocols

There are two possible approaches to a Gaussian semantics of protocols, that is, to a semantics characterized by keeping track of mean and variance of sample concentrations. They differ in the interpretation of the source of noise that produces the variances. We discuss the two options, and then choose one of the two based on both semantic simplicity and relevance to the well-mixed-solution assumption of chemical kinetics.

The first approach we call extrinsic noise. The protocol operations are deterministic, and the the evolution of each sample is also deterministic according to the rate equation (that is, Definition 4 (1)). Here, we imagine running the protocol multiple times over a distribution of initial conditions (i.e., the noise is given extrinsically to the evolution). The outcome is a final distribution that is determined uniquely by the initial distribution and by the deterministic evolution of each run. For example, the sample-split operation here is interpreted as follows. In each run, we have a sample whose concentration in general differs from its mean concentration over all runs. The two parts of a split are assigned identical concentration, which again may differ from the mean concentration. Hence, over all runs, the two parts of the split have identical variance, and their distributions are perfectly correlated, having the same concentration in each run. The subsequent time evolution of the two parts of the split is deterministic and hence identical (for, say, a $50\%$ split). In summary, in the extrinsic noise approach, the variance on the trajectories is due to the distribution of deterministic trajectories for the different initial conditions.

The second approach we call intrinsic noise. The protocol operations are again deterministic, but the evolution in each separate sample is driven by the Chemical Master Equation (i.e., the noise is intrinsic to the evolution). Here, we imagine running the protocol many times on the same initial conditions, and the variance of each sample is due ultimately to random thermal noise in each run. This model assumes a Markovian evolution of the underlying stochastic system, implying that no two events may happen at the same time (even in separate samples). Furthermore, as usual, we assume that chemical solutions are well-mixed: the probability of a reaction (a molecular collision) is independent of the initial position of the molecules in the solution. Moreover, the well-mixture of solutions is driven by uncorrelated thermal noise in separate samples. Given two initially identical but separate samples, the first chemical reaction in each sample (the first “fruitful” collision that follows a usually large number of mixing collisions) is determined by uncorrelated random processes, and their first reaction cannot happen at exactly the same time. Therefore, in contrast to the extrinsic noise approach, a sample-split operation results in two samples that are completely uncorrelated, at least on the time scale of the first chemical reaction in each sample. In summary, in the intrinsic noise approach, the variance on the trajectories is due to the distribution of stochastic trajectories for identical initial conditions.

In both approaches, the variance computations for mix and for split are the same: mix uses the squared-coefficient law to determine the variance of two mixed samples, while split does not change the variance of the sample being split. The only practical difference is that in the intrinsic noise interpretation we do not need to keep track of the correlation between different samples: we take it to be always zero in view of the well-mixed-solution assumption. As a consequence, each sample needs to keep track of its internal correlations represented as a local $\Sigma$ matrix of size $\left|\mathcal{A}\right|\times \left|\mathcal{A}\right|$, but not of its correlations with other samples, which would require a larger global matrix of potentially unbounded size. Therefore, the intrinsic noise interpretation results in a vast simplification of the semantics (Definition 6) that does not require keeping track of correlations of concentrations across separate samples. In the rest of the paper, we consider only the intrinsic noise semantics.

Definition 6

(Gaussian Semantics of a Protocol—Intrinsic Noise). The intrinsic-noise Gaussian semantics ${\left[\left[P\right]\right]}^{\rho }\in Prot\times Env\to S$ of a protocol $P\in Prot$ for a CRN $\mathcal{C}=(\mathcal{A},\mathcal{R})$, under environment $\rho \in Env=(Var\cup Par)\to S$, for a fixed horizon H with no ill-posed time evolutions, denotes the final CRN state $(\mu ,\Sigma ,V,T)\in S$ (Definition 2) of the protocol and is defined inductively as follows:   

$$\begin{array}{cc}\hfill & {\left[\left[x\right]\right]}^{\rho }=\rho \left(x\right)\hfill \\ \hfill & {\left[\left[(p_1\dots p_{\left|\mathcal{A}\right|},r_V,r_T)\right]\right]}^{\rho }=({\left[\left[p_1\right]\right]}^{\rho }\dots {\left[\left[p_{\left|\mathcal{A}\right|}\right]\right]}^{\rho },0^{\left|\mathcal{A}\right|\times \left|\mathcal{A}\right|},r_V,r_T)\hfill \\ \hfill & {\left[[letx=P_{1}in{P}_2]\right]}^{\rho }={\left[\left[P_2\right]\right]}^{{\rho }_1}\hfill \\ \hfill & where{\rho }_1=\rho \{x\leftarrow {\left[\left[P_1\right]\right]}^{\rho }\}\hfill \\ \hfill & {\left[\left[Mix(P_1,P_2)\right]\right]}^{\rho }=(\frac{V_1{\mu }_1+V_2{\mu }_2}{V_1+V_2},\frac{V_1^2{\Sigma }_1+V_2^2{\Sigma }_2}{{(V_1+V_2)}^2},V_1+V_2,\frac{V_1{T}_1+V_2{T}_2}{V_1+V_2})\hfill \\ \hfill & where({\mu }_1,{\Sigma }_1,V_1,T_1)={\left[\left[P_1\right]\right]}^{\rho }and({\mu }_2,{\Sigma }_2,V_2,T_2)={\left[\left[P_2\right]\right]}^{\rho }\hfill \\ \hfill & {\left[[letx,y=Split(P_1,p)in{P}_2]\right]}^{\rho }={\left[\left[P_2\right]\right]}^{{\rho }_1}\hfill \\ \hfill & wherer={\left[\left[p\right]\right]}^{\rho },0<r<1and(\mu ,\Sigma ,V,T)={\left[\left[P_1\right]\right]}^{\rho }\hfill \\ \hfill & and{\rho }_1=\rho \{x\leftarrow (\mu ,\Sigma ,rV,T),y\leftarrow (\mu ,\Sigma ,(1-r)V,T)\}\hfill \\ \hfill & {\left[\left[Equilibrate(P,p)\right]\right]}^{\rho }=({\mathbf{\mu }}_{\mu }\left(t\right),{\Sigma }_{\mu ,\Sigma }\left(t\right),V,T)\hfill \\ \hfill & wheret={\left[\left[p\right]\right]}^{\rho }and(\mu ,\Sigma ,V,T)={\left[\left[P\right]\right]}^{\rho }\hfill \\ \hfill & {\left[\left[Dispose\left(P\right)\right]\right]}^{\rho }=(0^{\left|\mathcal{A}\right|},0^{\left|\mathcal{A}\right|\times \left|\mathcal{A}\right|},0,0)\hfill \end{array}$$

together with ${\left[\left[p\right]\right]}^{\rho }$ defined as:

$$\begin{array}{cc}\hfill & {\left[\left[z\right]\right]}^{\rho }=\rho \left(z\right)\hfill \\ \hfill & {\left[\left[r\right]\right]}^{\rho }=r\hfill \end{array}$$

The semantics of $Equilibrate$ derives from Definition 4. The substitution notation $\rho \{x\leftarrow v\}$ represents a function that is identical to ρ except that at x it yields v; this may be extended to the case where x is a vector of distinct variables and v is a vector of equal size. The variables introduced by $let$ are used linearly (i.e., must occur exactly once in their scope), implying that each sample is consumed whenever used.

We say that the components $(\mu ,\Sigma )\in {\mathbb{R}}_{\ge 0}^{\left|\mathcal{A}\right|}\times {\mathbb{R}}^{\left|\mathcal{A}\right|\times \left|\mathcal{A}\right|}$ of a CRN state $(\mu ,\Sigma ,V,T)\in S$ form a Gaussian state, and sometimes we say that the protocol semantics produces a Gaussian state (as part of a CRN state).

The semantics of Definition 6 combines the end states of sub-protocols by linear operators, as we show in the examples below. We stress that it does not, however, provide a time-domain GP: just the protocol end state as a Gaussian state (random variable). In particular, a protocol like $letx=Equilibrate(P,p)inDispose\left(x\right)$ introduces a discontinuity at the end time p, where all the concentrations go to zero; other discontinuities can be introduced by $Mix$. A protocol may have a finite number of such discontinuities corresponding to liquid handling operations, and otherwise proceeds by LNA-derived GPs and by linear combinations of Gaussian states at the discontinuity points.

It is instructive to interpret our protocol operators as linear operators on Gaussian states: this justifies how covariance matrices are handled in Definition 6, and it easily leads to a generalized class of possible protocol operators. We discuss linear protocol operators in the examples below.

Example 2.

The $Mix$ operator combines the concentrations, covariances, and temperatures of the two input samples proportionally to their volumes. The resulting covariance matrix, in particular, can be justified as follows. Consider two (input) states $A=({\mu }_A,{\Sigma }_A,V_A,T_A),B=({\mu }_B,{\Sigma }_B,V_B,T_B)$ with $|{\mu }_A|=|{\mu }_B|=k$ and $|{\Sigma }_A|=|{\Sigma }_B|=k\times k$. Let 0 and 1 denote null and identity vectors and matrices of size k and $k\times k$. Consider a third null state $C=(0,0,V_C,T_C)$, with $V_C=V_A+V_B$ and $T_C=\frac{V_A{T}_A+V_B{T}_B}{V_C}$. The joint Gaussian distribution of $A,B,C$ is given by $\mu =\left[\begin{array}{c}{\mu }_A\\ {\mu }_B\\ 0\end{array}\right],$ $\Sigma =\left[\begin{array}{ccc}{\Sigma }_A& 0& 0\\ 0& {\Sigma }_B& 0\\ 0& 0& 0\end{array}\right].$ Define the symmetric hollow linear operator $Mix=\left[\begin{array}{ccc}0& 0& \frac{V_A}{V_C}1\\ 0& 0& \frac{V_B}{V_C}1\\ \frac{V_A}{V_C}1& \frac{V_B}{V_C}1& 0\end{array}\right].$ The zeros on the diagonal (hollowness) imply that the inputs states are zeroed after the operation, and hence discarded. Applying this operator to the joint distribution, by the linear combination of normal random variables we obtain a new Gaussian distribution with:

$$\begin{array}{cc}\hfill & {\mu }^{\prime }=Mix·\mu =\left[\begin{array}{c}0\\ 0\\ \frac{V_A{\mu }_A+V_B{\mu }_B}{V_C}\end{array}\right]\hfill \\ \hfill & {\Sigma }^{\prime }=Mix·\Sigma ·Mi{x}^{\top }=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& \frac{V_A^2{\Sigma }_A+V_B^2{\Sigma }_B}{V_C^2}\end{array}\right]\hfill \end{array}$$

Hence, all is left of ${\mu }^{\prime },{\Sigma }^{\prime }$ is the output state $C=({\mu }_C,{\Sigma }_C,V_C,T_C)$ where ${\mu }_C=\frac{V_A{\mu }_A+V_B{\mu }_B}{V_C}$ and ${\Sigma }_C=\frac{V_A^2{\Sigma }_A+V_B^2{\Sigma }_B}{V_C^2}$ as in Definition 6.

Example 3.

The $Split$ operator splits a sample in two parts, preserving the concentration and consequently the covariance of the input sample. The resulting covariance matrix, in particular, can be justified as follows. Consider one (input) state $A=({\mu }_A,{\Sigma }_A,V_A,T_A)$ and two null states $B=(0,0,p{V}_A,T_A)$, $C=(0,0,(1-p)V_A,T_A)$ with $0<p<1$. As above, consider the joint distribution of these three states, $\mu =\left[\begin{array}{c}{\mu }_A\\ 0\\ 0\end{array}\right],$ $\Sigma =\left[\begin{array}{ccc}{\Sigma }_A& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ and define the symmetric hollow linear operator $Split=\left[\begin{array}{ccc}0& 1& 1\\ 1& 0& 0\\ 1& 0& 0\end{array}\right].$ The 1s imply that concentrations are not affected. The whole submatrix of output states being zero implies that any initial values of output states are ignored. Applying this operator to the joint distribution, we obtain:

$$\begin{array}{cc}\hfill & {\mu }^{\prime }=Split·\mu =\left[\begin{array}{c}0\\ {\mu }_A\\ {\mu }_A\end{array}\right]\hfill \\ \hfill & {\Sigma }^{\prime }=Split·\Sigma ·Spli{t}^{\top }=\left[\begin{array}{ccc}0& 0& 0\\ 0& {\Sigma }_A& {\Sigma }_A\\ 0& {\Sigma }_A& {\Sigma }_A\end{array}\right]\hfill \end{array}$$

By projecting ${\mu }^{\prime },{\Sigma }^{\prime }$ on B and C, we are left with the two output states $B=({\mu }_A,{\Sigma }_A,p{V}_A,T_A)$, $C=({\mu }_A,{\Sigma }_A,(1-p)V_A,T_A)$, as in Definition 6. The correlation between B and C that is present in ${\Sigma }^{\prime }$ is not reflected in these outcomes: it is discarded on the basis of the well-mixed-solution assumption.

Example 4.

In general, any symmetric hollow linear operator, where the submatrix of the intended output states is also zero, describes a possible protocol operation over samples. As a further example, consider a different splitting operator, $letx,y=Osmo(P_1,p)in{P}_2$, that intuitively works as follows. A membrane permeable only to water is placed in the middle of an input sample $A=({\mu }_A,{\Sigma }_A,V_A,T_A)$, initially producing two samples of volume $V_A/2$, but with unchanged concentrations. Then, an osmotic pressure is introduced (by some mechanism) that causes a proportion $0<p<1$ of the water (and volume) to move from one side to the other, but without transferring any other molecules. As the volumes change, the concentrations increase in one part and decrease in the other. Consider, as in $Split$, two initial null states B and C and the joint distributions of those three states $\mu ,\Sigma$. Define a symmetric linear operator $Osmo=\left[\begin{array}{ccc}0& 1-p& p\\ 1-p& 0& 0\\ p& 0& 0\end{array}\right].$ Applying this operator to the joint distribution we obtain:

$$\begin{array}{cc}\hfill & {\mu }^{\prime }=Osmo·\mu =\left[\begin{array}{c}0\\ (1-p){\mu }_A\\ p{\mu }_A\end{array}\right]\hfill \\ \hfill & {\Sigma }^{\prime }=Osmo·\Sigma ·Osm{o}^{\top }=\left[\begin{array}{ccc}0& 0& 0\\ 0& {(1-p)}^2{\Sigma }_A& p(1-p){\Sigma }_A\\ 0& p(1-p){\Sigma }_A& p^2{\Sigma }_A\end{array}\right]\hfill \end{array}$$

producing the two output states $B=((1-p){\mu }_A,{(1-p)}^2{\Sigma }_A,p{V}_A,T_A)$ and $C=(p{\mu }_A,$$p^2{\Sigma }_A,(1-p)V_A,T_A)$ describing the situation after the osmotic pressure is applied. Again, we discard the correlation between B and C that is present in ${\Sigma }^{\prime }$ on the basis of the well-mixed-solution assumption.

Example 5.

We compute the validity of some simple equivalences between protocols by applying the Gaussian Semantics to both sides of the equivalence. Consider the CRN state $Poisson(k,V,T)=\left(\right[k],[\left[k\right]],V,T)$ over a single species having mean $\left[k\right]$ (a singleton vector) and variance $\left[\right[k\left]\right]$ (a $1\times 1$ matrix). This CRN state can be produced by a specific CRN [19], but here we add it as a new primitive protocol and extend our semantics to include ${\left[\left[Poisson(k,V,T)\right]\right]}^{\rho }=(\left[k\right],\left[\left[k\right]\right],V,T)$. Note first that mixing two uncorrelated Poisson states does not yield a Poisson state (the variance changes to $\left[\right[k/2\left]\right]$). Then, the following equation holds:

$$\begin{array}{cc}\hfill & {\left[\left[Mix(Poisson(k,V,T),Poisson(k,V,T))\right]\right]}^{\rho }\hfill \\ \hfill =& {\left[[letx,y=Split(Poisson(k,V,T),0.5)inMix(x,y)]\right]}^{\rho }\hfill \\ \hfill =& \left(\right[k],[[k/2]],V,T)\hfill \end{array}$$

That is, in the intrinsic-noise semantics, mixing two correlated (by Split) Poisson states is the same as mixing two uncorrelated Poisson states, because of the implicit decorrelation that happens at each $Split$. An alternative semantics that would take into account a correlation due to $Split$ could instead possibly satisfy the appealing equation ${\left[[letx,y=Split(P,0.5)inMix(x,y)]\right]}^{\rho }={\left[\left[P\right]\right]}^{\rho }$ for any P (i.e., splitting and remixing makes no change), but this does not hold in our semantics for $P=Poisson(k,V,T)$, where the left hand side yields $\left(\right[k],[[k/2]],V,T)$ (i.e., splitting, decorrelating, and remixing makes a change).

The $Mix$ operator can be understood as diluting the two input samples into the resulting larger output volume. We can similarly consider a $Dilute$ operator that operates on a single sample, and we can even generalize it to concentrate a sample into a smaller volume. In either case, let W be a new forced volume, and U be a new forced temperature for the sample. Let us define ${\left[\left[Dilute(P,W,U)\right]\right]}^{\rho }=(\frac{V}{W}\mu ,\frac{V^2}{W^2}\Sigma ,W,U)$ where $(\mu ,\Sigma ,V,T)={\left[\left[P\right]\right]}^{\rho }$. Then, the following equation holds:

$$\begin{array}{cc}\hfill & {\left[\left[Mix(Dilute(P_1,W_1,U_1),Dilute(P_2,W_2,U_2))\right]\right]}^{\rho }\hfill \\ \hfill =& {\left[\left[Dilute(Mix(P_1,P_2),W_1+W_2,\frac{W_1{U}_1+W_2{U}_2}{W_1+W_2})\right]\right]}^{\rho }\hfill \\ \hfill =& (\frac{V_1{\mu }_1+V_2{\mu }_2}{W_1+W_2},\frac{V_1^2{\Sigma }_1+V_2^2{\Sigma }_2}{{(W_1+W_2)}^2},W_1+W_2,\frac{W_1{U}_1+W_2{U}_2}{W_1+W_2})\hfill \end{array}$$

Note that, by diluting to a fixed new volume W, we obtain a protocol equivalence that does not mention, in the equivalence itself, the volumes $V_i$ resulting from the sub-protocols $P_i$ (which of course occur in the semantics). If instead we diluted by a factor (e.g., factor = 2 to multiply volume by 2), then it would seem necessary to refer to the $V_i$ in the equivalence.

To end this section, we provide an extension of syntax and semantics for the optimization of protocols. An optimize protocol directive identifies a vector of variables within a protocol to be varied in order to minimize an arbitrary cost function that is a function of those variables and of a collection of initial values. The semantics can be given concisely but implicitly in terms of an $argmin$ function: a specific realization of this construct is then the subject of the next section.

Definition 7

(Gaussian Semantics of a Protocol—Optimization). Let $z\in Pa{r}^N$ be a vector of (optimization) variables, and $k\in Pa{r}^M$ be a vector of (initialization) variables with (initial values) $r\in {\mathbb{R}}^M$, where all the variables are distinct and $Par$ is disjoint from $Var$. Let $P\in Prot$ be a protocol for a CRN $\mathcal{C}=(\mathcal{A},\mathcal{R})$, with free variables z and k and with any other free variables covered by an environment $\rho \in Env=(Var\cup Par)\to S$. Given a cost function $C:{\mathbb{R}}^{\left|\mathcal{A}\right|}\times {\mathbb{R}}^N\to \mathbb{R}$, and a dataset $\mathcal{D}{\subseteq }_{fin}{\mathbb{R}}^M\times {\mathbb{R}}^N\times {\mathbb{R}}^{\left|\mathcal{A}\right|}$, the instruction “optimize k=r, z in P with C given $\mathcal{D}$" provides a vector in ${\mathbb{R}}^N$ of optimized values for the z variables (For simplicity, we omit extending the syntax with new representations for C and $\mathcal{D}$, and we let them represent themselves). The optimization is based on a stochastic process $\mathbf{X}\in {\mathbb{R}}^{M+N}\to {\mathbb{R}}_{\ge 0}^{\left|\mathcal{A}\right|}\times {\mathbb{R}}^{\left|\mathcal{A}\right|\times \left|\mathcal{A}\right|}$ derived from P via Definition 6.

$$\begin{array}{cc}\hfill & {\left[[optimizek=r,zinPwithCgiven\mathcal{D}]\right]}^{\rho }=\hfill \\ \hfill & argmi{n}_{u\in {\mathbb{R}}^N}{\mathbb{E}}_{y\sim \mathbf{X}(r,u|\mathcal{D})}\left[C(y,u)\right]\hfill \\ \hfill & where\mathbf{X}(r,u)={\left[\left[P\right]\right]}^{\rho \{k\leftarrow r,z\leftarrow u\}}forallr\in {\mathbb{R}}^{M}andu\in {\mathbb{R}}^N,\hfill \end{array}$$

where ${\mathbb{E}}_{y\sim \mathbf{X}(r,u|\mathcal{D})}[·]$ stands for the expectation with respect the conditional distribution of $\mathbf{X}(r,u)$ given $\mathcal{D}$.

From Definition 6, it follows that, for any $(r,u)\in {\mathbb{R}}^{M+N},$ $\mathbf{X}(r,u)$ is a Gaussian random variable. In what follows, for the purpose of optimization, we further assume that $\mathbf{X}$ is a Gaussian process defined on the sample space ${\mathbb{R}}^{M+N}$, i.e., that for any possible $(r_1,u_1)$ and $(r_2,u_2),$$\mathbf{X}(r_1,u_1)$ and $\mathbf{X}(r_2,u_2)$ are jointly Gaussian. Such an assumption is standard [6] and natural for our setting (e.g., it is guaranteed to hold for different equilibration times) and guarantees that we can optimize protocols by employing results from Gaussian process optimization as described in detail in Section 2.4.

2.4. Optimization of Protocols through Gaussian Process Regression

In Definition 7 we introduced the syntax and semantics for the optimization of a protocol from data. In this section, we show how the optimization variables can be selected in practice. In particular, we leverage existing results for Gaussian processes (GPs) [6]. We start by considering the dataset $\mathcal{D}=\{(r_i,u_i,y_i),i\in \{1,\dots ,n_{\mathcal{D}}\}\}{\subseteq }_{fin}{\mathbb{R}}^M\times {\mathbb{R}}^N\times {\mathbb{R}}^{\left|\mathcal{A}\right|}$ comprising $n_{\mathcal{D}}$ executions of the protocol for possibly different initial conditions, i.e., each entry gives the concentration of the species after the execution of the protocol $\left(y_i\right)$, where the protocol has been run with $k=r_i,z=u_i$. Then, we can predict the output of the protocol starting from $\overline{x}=(r,u)$ by computing the conditional distribution of $\mathbf{X}$ (Gaussian process introduced in Definition 7) given the data in $\mathcal{D}$. In particular, under the assumption that $\mathbf{X}$ is Gaussian, it is well known that the resulting posterior model, $\mathbf{X}\left(\overline{x}\right|\mathcal{D})$, is still Gaussian with mean and covariance functions given by

$$\begin{array}{cc}\hfill & {\mathbf{\mu }}_p\left(\overline{x}\right)=\mathbf{\mu }\left(\overline{x}\right)+{\Sigma }_{\overline{x},\mathcal{D}}{({\Sigma }_{\mathcal{D},\mathcal{D}}+{\sigma }^{2}I)}^{-1}(y_{\mathcal{D}}-{\mathbf{\mu }}_{\mathcal{D}})\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill & {\Sigma }_p(\overline{x},\overline{x})={\Sigma }_{\overline{x},\overline{x}}-{\Sigma }_{\overline{x},\mathcal{D}}{({\Sigma }_{\mathcal{D},\mathcal{D}}+{\sigma }^{2}I)}^{-1}{\Sigma }_{\overline{x},\mathcal{D}}^T,\hfill \end{array}$$

(4)

where ${\sigma }^{2}I$ is a diagonal covariance modeling i.i.d. Gaussian observation noise with variance ${\sigma }^2$, $\mathbf{\mu }\left(\overline{x}\right)$ and ${\Sigma }_{\overline{x},\overline{x}}$ are the prior mean and covariance functions, ${\Sigma }_{\mathcal{D},\overline{x}}$ is the covariance between $\overline{x}$ and all the points in $\mathcal{D}$, and $y_{\mathcal{D}}$, ${\mathbf{\mu }}_{\mathcal{D}}$ are vectors of dimensions ${\mathbb{R}}^{\left|\Lambda \right|·n_{\mathcal{D}}}$ containing for all $(x_i,u_i,y_i)\in \mathcal{D}$, respectively, $y_i$ and $\mathbf{\mu }\left((x_i,u_i)\right)$ [6]. Note that, in order to encode our prior information of the protocol, we take $\mathbf{\mu }\left(\overline{x}\right)$ to be the mean of the GP as defined in Definition 6, while for the variance we can have ${\Sigma }_{\overline{x},\overline{x}}$ to be any standard kernel [6]; in the experiments, as is standard in the literature, we consider the widely used squared exponential kernel, which is expressive enough to approximate any continuous function arbitrarily well [20]. However, we remark that, for any parametric kernel, such as the squared exponential, we can still select the hyper-parameters that best fit the variance given by Definition 6 as well as the data [6]. We should also stress that the resulting $\mathbf{X}\left(\overline{x}\right|\mathcal{D})$ is a Gaussian process that merges in a principled way (i.e., via Bayes’ rule) our prior knowledge of the model (given by Definition 6) with the new information given in $\mathcal{D}.$ Furthermore, it is worth stressing again that $\mathbf{X}\left(\overline{x}\right|\mathcal{D})$ is a GP with input space given by ${\mathbb{R}}^{M+N},$ which is substantially different from the GP defined by the LNA (Definition 5), where the GP was defined over time.

For a given set of initialization parameters r our goal is to synthesize the optimization variables that optimize the protocol with respect to a given cost specification $C:{\mathbb{R}}^{\left|\mathcal{A}\right|}\times {\mathbb{R}}^N\to \mathbb{R}$, i.e., we want to find $u^{*}$ such that

$$\begin{array}{c}\hfill u^{*}=argmi{n}_{u\in {\mathbb{R}}^N}{\mathbb{E}}_{y\sim \mathcal{N}({\mathbf{\mu }}_p\left(\overline{x}\right),{\Sigma }_p(\overline{x},\overline{x}))}\left[C(y,u)\right],\end{array}$$

(5)

where ${\mathbb{E}}_{y\sim \mathcal{N}({\mathbf{\mu }}_p\left(\overline{x}\right),{\Sigma }_p(\overline{x},\overline{x}))}[·]$ is the expectation with respect to the GP given by Equations (3) and (4)). Note that r is a known vector of reals (see Definition 7), hence we only need to optimize for the free parameters u. In general, the computation of $u^{*}$ in Equation (5) requires solving a non-convex optimization problem that cannot be solved exactly [21]. In this paper, we approximate Equation (5) via gradient-based methods. These methods, such as gradient descent, require the computation of the gradient of the expectation of C with respect to u:

$$\begin{array}{c}\hfill \frac{\partial {\mathbb{E}}_{y\sim \mathcal{N}({\mathbf{\mu }}_p\left(\overline{x}\right),{\Sigma }_p(\overline{x},\overline{x}))}\left[C(y,u)\right]}{\partial u}.\end{array}$$

(6)

Unfortunately, direct computation of the gradient in (6) is infeasible in general, as the probability distribution where the expectation is taken depends on u itself (note that $\overline{x}=(r,u)$). However, for the GP case, as shown in Lemma 1, the gradient of interest can be computed directly by reparametrizing the Gaussian distribution induced by Equations (3) and (4).

Lemma 1.

For $\overline{x}=(r,u)$, let $D\left(\overline{x}\right)$ be the matrix such that $D\left(\overline{x}\right)D^T\left(\overline{x}\right)={\Sigma }_p(\overline{x},\overline{x})$. Then, it holds that

$$\frac{\partial {\mathbb{E}}_{y\sim \mathcal{N}({\mathbf{\mu }}_p\left(\overline{x}\right),{\Sigma }_p(\overline{x},\overline{x}))}\left[C(y,u)\right]}{\partial u}={\mathbb{E}}_{z\sim \mathcal{N}(0,I)}\left[\frac{\partial C({\mathbf{\mu }}_p\left(\overline{x}\right)+D\left(\overline{x}\right)z,u)}{\partial u}\right].$$

(7)

In particular, $D\left(\overline{x}\right)$ as defined above is guaranteed to exist under the assumption that ${\Sigma }_p$ is positive definite and can be computed via Cholesky decomposition. Furthermore, note that if the output is uni-dimensional, then $D\left(\overline{x}\right)$ is simply the standard deviation of the posterior GP in $\overline{x}$.

Example 6.

Given the CRN $\mathcal{C}$ introduced in Example 1, consider the protocol

$$P=Equilibrate\left(\right(0.1mM,0.001mM,0.001mM,1\mu L,20C),T),$$

which seeks to evolve the various species from an initial starting concentration of $0.1$ mM for species a and $0.001$ mM for both species b and c for T seconds. Hence, for this example, we have $r=[0.1,0.001,0.001,1,20]$ including initial conditions, volume, and temperature and $u=\left[T\right]$, that is, the only optimization variable is the equilibration time. We consider a dataset $\mathcal{D}$ composed of the following six data points (for simplicity we report the output values of only species b and omit unit of measurement, which are mM for concentration, μL for volume, Celsius degrees for temperature, and seconds for equilibration time):

$$\begin{array}{cc}\hfill & d_1=((0.1001,0.0015,0.001,1,20),0,0.001)d_2=((0.099,0.001,0.001,1,20),40,0.001)\hfill \\ \hfill & d_3=((0.1,0.001,0.001,1,20),150,0.09)d_4=((0.1,0.001,0.002,1,20),250,0.08)\hfill \\ \hfill & d_5=((0.09,0.001,0.0015,1,20),400,0.003)d_6=((0.1,0.001,0.002,1,20),500,0.001),\hfill \end{array}$$

where, for example, in $d_1$ we have that the initial concentration for $a,b,c$ are, respectively, $0.1001$, $0.0015$, $0.001$, volume and temperature at which the experiment is performed are 1 and 20 and the observed value for b at the end of the protocol for $T=0$ is $0.001.$ We assume an additive Gaussian observation noise (noise in the collection of the data) with standard deviation $\sigma =0.01$. Note that the initial conditions of the protocol in the various data points in general differ from those of P, which motivates having Equations (3) and (4) dependent on both r and u.

We consider the prior mean given by Equation (1) and independent squared exponential kernels for each species with hyper-parameters optimized through maximum likelihood estimation (MLE) [6]. The resulting posterior GP is reported in Figure 3, where is compared with the prior mean and the true underlying dynamics for species b (assumed for this example to be deterministic). We note that with just 6 data points in $\mathcal{D}$, the posterior GP is able to correctly recover the true behavior of b with relatively low uncertainty.

We would like to maximize the concentration of b after the execution of the protocol. This can be easily encoded with the following cost function:

$$C(a,b,c,u)=-b^2.$$

For $r=[0.1,0.001,0.001,1,20]$, the value obtained is $T\sim 230$, which is simply the one that maximizes the posterior mean. In Section 3.1, we will show how, for more complex specifications, the optimization process will become more challenging because of the need to also account for the variance in order to balance between exploration and exploitation.

3. Results

We consider two case studies where we illustrate the usefulness of our framework. The first example illustrates how our framework can be employed to optimize a Gibson assembly protocol [9] from experimental data. In the second example, we illustrate the flexibility of our language and semantics on a combination of protocol operations. Furthermore, we also use this example to show how our framework can be employed to perform analysis of the protocol parameters while also accounting for the uncertainty in both the model dynamics, the protocol parameters, and the data collection.

3.1. Gibson Assembly Protocol

We start by considering a simplified version of the Gibson assembly, a protocol widely used for joining multiple DNA fragments [9]. We consider the following model of Gibson assembly described by Michaelis–Menten kinetics, where two DNA double strands, AB and BA (with shared endings A and B, but different middle parts), are enzymatically joined in two symmetric ways according to their common endings, and the resulting strands ABA and BAB are circularized into the same final structure O. The resulting dynamical model is given by the following ODEs (assuming AB is over abundant with respect to BA):

$$\begin{array}{cc}\hfill & \frac{dAB\left(t\right)}{dt}=0\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & \frac{dBA\left(t\right)}{dt}=-\frac{k_{ca{t}_1}·AB\left(t\right)·BA\left(t\right)}{BA\left(t\right)+K_{m_1}}-\frac{k_{ca{t}_2}·AB\left(t\right)·BA\left(t\right)}{BA\left(t\right)+K_{m_2}}\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & \frac{dABA\left(t\right)}{dt}=\frac{k_{ca{t}_1}·AB\left(t\right)·BA\left(t\right)}{BA\left(t\right)+K_{m_1}}-\frac{k_{ca{t}_1}·ABA\left(t\right)}{ABA\left(t\right)+K_{m_1}}\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & \frac{dBAB\left(t\right)}{dt}=\frac{k_{ca{t}_2}·AB\left(t\right)·BA\left(t\right)}{BA\left(t\right)+K_{m_2}}-\frac{k_{ca{t}_2}·BAB\left(t\right)}{BAB\left(t\right)+K_{m_2}}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & \frac{dO\left(t\right)}{dt}=\frac{k_{ca{t}_1}·ABA\left(t\right)}{ABA\left(t\right)+K_{m_1}}+\frac{k_{ca{t}_2}·BAB\left(t\right)}{BAB\left(t\right)+K_{m_2}},\hfill \end{array}$$

(12)

where $k_{ca{t}_1},k_{ca{t}_2},K_{m_1},K_{m_2}$ are given parameters. Gibson assembly is a single-step isothermal protocol: once all the ingredients are combined, it can be written as:

$$P=Equilibrate(([1mM,x_{BA}mM,0mM,0mM,0mM],1\mu L,50C),T),$$

where $[1,x_{BA},0,0,0]$ is a vector of the initial concentration of the various species with $AB$ initialized at 1 mM, $BA$ at $x_{BA}$ mM, and all the other species to 0 mM for a sample of volume of 1 $\mathsf{\mu }$L and at a temperature 50 Celsius degrees, equilibrated for T seconds. For this example, we have that the set of optimization variables is $u=[x_{BA},T]$ and the goal is to find initial condition of $BA$ ($x_{BA}$) and equilibration time (T) such that, at the end of the protocol, species O has a concentration as close as possible to a desired value $O_{des},$ while also keeping the equilibration time and the initial concentration of $BA$ close to some reference values. This can be formalized with the following cost function:

$$C(x_{BA},T)={\left(\beta (O\left(T\right)-O_{des})\right)}^2+\lambda {\left(\frac{T-T_{ref}}{T_{norm}}\right)}^2+(1-\lambda ){(x_{BA}-I_{ref})}^2,$$

where $O\left(T\right)$ is the concentration of strand O after the execution of the protocol. $T_{ref},$$T_{norm},$$I_{ref},$ and $\beta ,$ are parameters describing, respectively, the reference equilibration time, a normalization term for the equilibration time, the reference initial concentration of $BA$, and a weight term such that if $\beta >1$ we give more importance to $O\left(T\right)$ being close to $O_{des}$ compared to equilibration time and initial conditions being closer to their reference values. Similarly, $\lambda \in [0,1]$ is a weight term that balances the importance that equilibration time and the initial concentration of $BA$ are close to their reference values. In our experiments, we fix $O_{des}=0.6,$$T_{ref}=700,$$T_{norm}=1000,$$I_{ref}=1$, and $\beta =20.$

The dynamical model considered in Equations (8)–(12) is obtained by combining multiple enzymatic reactions into single Michaelis–Menten steps and the value of the reaction rates may not be accurate. Hence, to obtain a more accurate model, we combine Equations (8)–(12) with the experimental data from [9] under the assumption that observations are corrupted by a Gaussian noise of relatively high standard deviation of $0.1.$

In Figure 4, we plot the synthesized (approximately) optimal values of T and $x_{BA}$ for various values of $\lambda$ by using gradient descent, with the gradient computed as shown in Lemma 1. Note that even for $\lambda =0$ or $\lambda =1$ the cost cannot be made identically zero. This is due to the uncertainty in the model that leads to a non-zero variance everywhere. Note also that for $\lambda =1$ the synthesized time T is smaller than 700 (the reference equilibration time). This can be explained by looking at Figure 4 (Right), where we plot the variance of $O\left(T\right)$ as a function of T and $x_{BA}$. In fact, as we have available only data (reported in the Appendix B) for $x_{BA}\sim 1$ the variance increases when $x_{BA}$ is far from 1. As a consequence, the algorithm automatically balances this tradeoff by picking an equilibration time that is close to the target value, but also allows for an $x_{BA}$ close to 1 in order the keep the variance under control.

3.2. Split and Mix Protocol

Consider the CRN $\mathcal{C}=\left(\right\{a,b,c\},\mathcal{R})$, where $\mathcal{R}$ is given by the reactions:

$$a+b{\to }^{1}b+ca+c{\to }^{1}a+ab+c{\to }^{1}c+c.$$

Consider also the associated protocol $P_{split\&mix}$ shown below in the syntax of Definition 5 (with underscore standing for an unused variable, with initial concentration vectors ordered as $(a_0,b_0,c_0)$, and initial states matching the pattern $((a_0,b_0,c_0)$, $volume$, $temperature)$:

$$\begin{array}{cc}\hfill & letA=\left(\right(10mM,0mM,1mM),1\mu L,20C)in\hfill \\ \hfill & letA1=Equilibrate(A,100s)in\hfill \\ \hfill & letC,D=Split(A1,0.5)in\hfill \\ \hfill & let\_=Dispose\left(C\right)in\hfill \\ \hfill & letB=\left(\right(0mM,10mM,1mM),1\mu L,20C)in\hfill \\ \hfill & letB1=Equilibrate(B,100s)in\hfill \\ \hfill & letE=Mix(D,B1)in\hfill \\ \hfill & Equilibrate(E,1000s)\hfill \end{array}$$

This protocol interleaves equilibration steps with mixing and splitting. During the first $Equilibrate$, only the second reaction is active in sample A, due to its initial conditions, yielding monotonic changes in concentrations. During the second $Equilibrate$, similarly, only the third reaction is active in sample B. During the third $Equilibrate$, after samples A and B have been manipulated and mixed, all three reactions are active in sample E, yielding an oscillation.

The semantics of Definition 6 can be used to unravel the behavior of $P_{split\&mix}$. In particular, integration of the CRN $\mathcal{C}$ in samples A and B yields states $S_A$, $S_B$ that, after some split and mix manipulations, determine the initial conditions for both mean and covariance for the integration of $\mathcal{C}$ in sample E.

We can numerically evaluate the protocol following the semantics: this is shown in Figure 5 (Left and Center), where the protocol evaluation results in three simulations that implement the three integration steps, where each simulation is fed the results of the previous simulations and protocol operations.

To show the usefulness of the integrated semantics, in Figure 5 (Right) we perform a global sensitivity analysis of the whole protocol with respect to the protocol parameters $t_1=100$, $t_2=100$, $t_3=1000$, the duration of the three $Equilibrate$, and $s_1=0.5$, the proportion of $Split$. We produce density plots of the means and variances of a,b,c at the end of the protocol when $t_1$,$t_2$,$t_3$,$s_1$ vary jointly randomly by up to 10%. This is thus a sensitivity analysis of the joint effects of three simulations connected by liquid handling steps. It is obtained by providing the whole parameterized protocol, not its separate parts, to the harness that varies the parameters and plots the results. We could similarly vary the initial concentrations and reaction rates, and those together with the protocol parameters.

4. Discussion

Automation already helps scaling up the production of chemical and biochemical substances, but it is also hoped it will solve general reproduceability problems. In the extreme, even one-off experiments that have no expectation of reproduction or scaling-up should be automated, so that the provenance of the data they produce can be properly recorded. Most pieces of lab equipment are already computer controlled, but automating their interconnection and integration is still a challenge that results in poor record keeping. A large amount of bottom-up development will be necessary to combine all these pieces of equipment into “fabrication lines”. However, it is also important to start thinking top-down at what the resulting connecting “glue” should look like, as we have attempted to do in this paper. This is because, if automation is the ultimate goal, then some pieces of equipment that are hard to automate should be discarded entirely.

Consider, for example, the Split operation from Definition 6, which separates a sample into two samples. That represents a vast array of laboratory procedures, from simple pipetting to microfluidic separation, each having its own tolerances and error distributions (which are well characterized and could be included into the language and the analysis). Despite all those variations, it is conceivable that a protocol compiler could take an abstract Split instruction, and a known piece of equipment, and automatically tailor the instruction for that equipment. Digital microfluidics is particularly appealing in this respect because of the relative ease and uniformity of such tailoring [22]. Therefore, one could decide to work entirely within digital microfluidics, perhaps helping to make that area more robust, and avoid other kinds of equipment.

Are there drawbacks to this approach? First, as we mentioned, equipment should be selected based on ease of automation: whole new lab procedures may need to be developed in replacement, along with standards for equipment automation. Second, protocol languages are more opaque than either sets of mathematical equations or sequences of laboratory steps. This needs to be compensated for with appropriate user interfaces for common use in laboratories.

About the second point, as a first step, we have separately produced an easily deployed app (Kaemika [15]) that supports the integrated language described here, although pragmatic details of the syntax diverge slightly. We have used it to run simulations and to produce Figure 2 and Figure 5 (see Appendix B). It uses a standard ODE solver for simulating the “Equilibrate” steps of the protocols, i.e., the chemical reaction kinetics. It uses the linear noise approximation (LNA) for stochastic simulations; that is, means and variances of concentrations are computed by an extended set of ODEs and solved with the same ODE solver along with the deterministic rate equations. The liquid handling steps of the protocols are handled following the semantics of Definition 6, noting that in the stochastic case means and variances are propagated into and out of successive equilibrate and liquid handling steps. The result is an interleaved simulation of equilibrate and liquid handing steps, which is presented as an interleaved graphical rendering of concentration trajectories and of droplet movements on a simulated microfluidic device. The sequence of protocol steps can be extracted as a graph, omitting the detailed kinetics. The kinetic equations operating at each step of the protocol can be extracted as well.

In this context, we should also stress that our proposed Gaussian semantics makes various assumptions: (1) we assume that molecular processes described by CRNs and liquid handling operations can be modeled by a Gaussian process, and (2) we assume that the collected data are corrupted by additive Gaussian noise. The former assumption is justified by the precision of lab equipment to perform liquid handling operations, whose uncertainty is generally very small, and by the Central Limit Theorem (CLT). The CLT guarantees that the solution of the Chemical Master Equation will converge to a Gaussian process in the limit of high number of molecules [23], as is common in wet lab experiments. In particular, as illustrated in Experiment 5.1 in [8], already a number of molecules of the order of hundreds for each species generally guarantees that a Gaussian approximation is accurate. It is obvious that, in case of experiments with single molecules, such an approximation may be inaccurate and a full treatment of the CME would be more appropriate [24]. The assumption that observations are corrupted by Gaussian additive noise is standard [25]. However, we acknowledge that in certain scenarios, non-Gaussian or multiplicative observation noise may be required. In this case, we would like to stress that Gaussian process regression can still be performed with success, at the price of introducing approximations on the computation of the posterior mean and variance [6].

As automation standards are developed for laboratory equipment, we will be able to target our language to such standards, and extend it or adapt it as needed. In this context, our future works include extension to both the syntax and semantics of our language to include common laboratory procedures that we have not investigated here, for example, changing the temperature of a sample (as in a thermal cycler) or its volume (as in evaporation or dilution). These are easy to add to our abstract framework, but each corresponds to a whole family of lab equipment that may need to be modeled and integrated in detail. Furthermore, we plan to extend the semantics to allow for more general stochastic processes that may be needed when modeling single molecules.

5. Conclusions

We have introduced a probabilistic framework that rigorously describes the joint handling of liquid manipulation steps and chemical kinetics, throughout the execution of an experimental protocol, with particular attention to the propagation of uncertainty and the optimization of the protocol parameters. A central contribution of this paper is the distinction between the intrinsic and extrinsic approach to noise, which leads to a much simplified semantics under a chemically justified assumption of well-mixed-solutions. The semantics are reduced to operations on Gaussian states, where any modeled or hypothetical protocol operation is a symmetric linear operator on Gaussian states.

The Gaussian process semantics approach is novel to this work, and is novel with respect to the Piecewise Deterministic Markov Process semantics in [5], which treats chemical evolution as deterministic. The semantics in this work are about a collection of deterministic protocol operations, but note that (1) stochasticity in chemical kinetics is propagated across protocol operations, making the whole protocol stochastic, and (2) we have shown examples of how to easily incorporate new protocol operators as linear operators on Gaussian states, which may include ones that introduce their own additional stochasticity. The Gaussian approach in this paper enables principled integration of a protocol model with experimental data, which in turn enables automated optimization and analysis of experimental biological protocols. The syntax of protocols is chosen for mathematical simplicity, reflecting the one in [5]; a closely related but more pragmatic syntax is now implemented in [15].