# Bayesian and Algebraic Strategies to Design in Synthetic Biology

Robyn P. Araujo, Sean T. Vittadello and Michael P.H. Stumpf R.P. Araujo is
with the Queensland University of TechnologyS. T. Vittadello and Michael P.H.
Stumpf are with the University of Melbourne

###### Abstract

Innovation in synthetic biology often still depends on large-scale
experimental trial-and-error, domain expertise, and ingenuity. The application
of rational design engineering methods promise to make this more efficient,
faster, cheaper and safer. But this requires mathematical models of cellular
systems. And for these models we then have to determine if they can meet our
intended target behaviour. Here we develop two complementary approaches that
allow us to determine whether a given molecular circuit, represented by a
mathematical model, is capable of fulfilling our design objectives. We discuss
algebraic methods that are capable of identifying general principles
guaranteeing desired behaviour; and we provide an overview over Bayesian
design approaches that allow us to choose from a set of models, that model
which has the highest probability of fulfilling our design objectives. We
discuss their uses in the context of biochemical adaptation, and then consider
how robustness can and should affect our design approach.

###### Index Terms:

Synthetic biology; approximate Bayesian computation; adaptation; robustness.

##  I Introduction

Trial and error coupled to ingenuity have arguably been responsible for many
important innovations in biotechnology. Unfortunately, however, these do not
scale. Instead recent progress in synthetic biology is generally based on
rational design engineering approaches [1, 2]. Design in cellular and
molecular biology cannot be based directly on the same physical and chemical
principles at work in technological systems. Naturally evolved systems are
just too complex; our knowledge of their structure and dynamics too primitive;
and our ability to interfere with and control their behaviour is still very
limited.

Here we discuss a set of complementary approaches that can be used in order to
design synthetic biological systems that exhibit certain types of behaviour,
or that fulfil certain types of biological functions. We need to acknowledge
the differences between biological and technical systems and our ignorance of
the details of the molecular mechanisms driving important cellular functions
and phenotypes. In light of this we need better ways of developing meaningful,
realistic, reliable, and validated mathematical models. And in the context of
synthetic biology we need quick and robust ways of picking those models which
can meet our design targets and objectives. The latter is the main purpose of
this paper.

Below we first outline the problem of design in synthetic biology: we are
interested in dynamical systems, implemented in synthetically engineered
molecular circuits, that fulfil certain design objectives. One overarching
problem is that they need to function robustly and reliably on top of the
cellular physiology. We use modelling or in silico analysis as a way of
triaging the vast potential design space in order to identify the best
designs, and we will attempt to outline this problem and the confounding
mathematical factors in the necessary detail, before illustrating two
complementary strategies of exploring this design space. One way of assessing
the promise of designs is to assess their level of robustness to uncertainties
in the models that we investigate, and we will make more precise what we mean
by robustness.

We will provide an overview of how algebraic approaches can be used to distill
design principles for, at least, certain types of behaviour. In particular, we
illustrate how such an approach helps to identify absolute requirements for a
dynamical system to exhibit a fundamentally important network response known
as Robust Perfect Adaptation (RPA) [3]. Algebraic approaches are capable of
yielding qualitative insights; Bayesian design methods provide an alternative
approach[4]. Unlike the algebraic approaches they are computationally costly,
but provide rational criteria for choosing from among competing models. We
illustrate in general terms how this can be used in practice, again in the
context of adaptive responses.

We conclude by a discussion how robustness of mathematical models can reassure
us that a design may function in practice. We argue that taken together,
algebraic methods, Bayesian design, and robustness, offer enormous scope for
enabling rational design in synthetic biology.

##  II Design in Synthetic Biology: Objectives and Strategies

In synthetic biology we are interested in identifying a (dynamical) system
that reliably produces a certain type of target behaviour. For example, for a
given input, $x$, our system state $y$ should produce a target output or
behaviour, $T$. Both $x$ and $T$, can be vector-valued though often they may
not be; but $y$, the system state, will be a vector, including all the
relevant molecular species (mRNA, proteins, post-translationally modified
proteins, metabolites, etc.) involved in a given process.

We consider systems of the following form,

$\frac{dy}{dt}=f(y;x,t,\theta)$ (1)

where $\theta$ denotes the system parameters (which are typically unknown). We
want to determine if the system given by $f(y;x,t,\theta)$ is capable of
producing the desired output, $T$, and meet the design objectives.

In our applications we may be interested in systems that can exhibit switch-
like behaviour, filter out high frequency noise, exhibit perfect adaptation,
or more complex phenotypes, such as the ability of a cellular system to
produce Turing pattern behaviour. Specifying the design objectives, $T$, can
involve qualitative characteristics, or quantitative aspects of system
behaviour. Often we may resort to surrogate data that captures the desired
behaviour; this is especially the case if precise mathematical or algebraic
definitions of system behaviour are difficult; or if there are no simple
qualitative features (such as algebraic identities).

Model-based design in synthetic biology can be done in different ways: (i) we
can use domain expertise to formulate mathematical models, $f(y;x,t,\theta)$,
determine if they fulfil our design objectives, and optimise the performance
by varying $\theta$ or by making small alterations to the model; (ii) choose
from a set of models $\mathcal{F}=\\{f_{1},f_{2},\ldots,f_{N}\\}$ using, for
example statistical model selection in order to determine which models are
best able to deliver the design objectives; (iii) explore exhaustively a class
of relevant models. That is, now $\mathcal{F}$ might contain all models, $f$,
with 2, 3, or 4 species or nodes (the dimension of the vector $y$) to
determine if they meet the design objectives; (iv) identify universal
characteristics guaranteeing the design objectives are met and construct
models (and, as the end-target, synthetic biological systems) embodying these
principles. The first approach relies a lot on luck in getting the initial
model right. The second and third approach can become computationally
prohibitively expensive even for small networks; they are, however, flexible
and widely applicable. The final approach is less ad hoc than the others, and
can provide guarantees on the behaviour of designs. But it requires a more
grounded theoretical understanding of the design objectives than we often
have.

Of course, there are also ways in which we can combine the last three
approaches to gain confidence in our analyses and designs, and to streamline
the design process. If, for example, design principles are known we can use
any such theoretical principles to reduce the set of potential candidate
models, $\mathcal{F}$, potentially quite considerably.

##  III Algebraic Approaches to Discovering Design Principles

In this section we consider algebraic approaches that may be able to identify,
within the vast potential design space, the exceedingly limited regions that
contain the ‘best’ designs – those that can perform an essential biological
function in a robust and reliable manner. We emphasize that such analytical
methods may not be applicable to the study of all possible design objectives,
and may not be well suited to the elucidation of complex biological behaviours
such as multistationarity [5] or Turing patterning[6]. But algebraic (and
ultimately topological) methods are beginning to show tremendous promise in
the context of robust asymptotic tracking problems, and have recently
delineated the full solution space for network topologies capable of Robust
Perfect Adaptation (RPA) [3]. RPA is is of fundamental significance in
biological systems, both natural and synthetic, and has been widely observed
in cellular signalling networks [7, 3, 8], sensory systems [9, 10], in the
regulation of bacterial chemotaxis [11, 12], bacterial stress response [13,
14], intracellular osmolarity [15], sigma-70 activity [16], and homeostatic
control of plasma mineral concentrations [17, 18], and is thought to play a
role in robust patterning during organism development [19, 20].

###  III-A The Control Theory viewpoint on RPA

RPA is the ability of a system to generate an output that returns to a fixed
reference level (its ‘setpoint’) following a persistent change in the system
input, with no need for parameter fine-tuning [3]. RPA is thus equivalent to
the well-established control theoretic notion of robust steady-state tracking
or constant disturbance rejection, while absolute concentration robustness
(ACR, a special type of RPA [21]) is equivalent to the control theoretic
concept of robustness to disturbances in initial conditions [22]. This
fundamental problem of robust asymptotic tracking of a desired trajectory,
while rejecting unwanted disturbances, is one of the defining problems of
classical automatic control [23]. In the 1970s, the landmark studies of
Francis and Wonham [24, 25] investigated the necessary controller structures
to achieve robust regulation of a system, with internal stability, and
established what is now known as the Internal Model Principle (IMP). By this
principle, a control system is able to reject exogenous stimuli or
disturbances by incorporating within itself a model of the dynamic structure
of the network stimuli. For RPA, then, the internal model must produce
constant signals and is equivalent to the requirement for integral (feedback)
control [26]. Thus, in all RPA-capable biochemical reaction networks, we
expect an integrator - an internal model of a system that produces constant
inputs or disturbances - to be embedded or concealed within the biochemical
reaction rates that govern the dynamical response of the network.

Suppose we have a $n$-dimensional system of interacting molecules, with state
vector $x(t)$ representing the concentrations of the molecules. In general, we
seek some change of coordinates that identifies a subsystem, $z$, of the form

$\frac{dz}{dt}=p(x)(\hat{x}-c),$ (2)

where $\hat{x}$ represents the concentration of the RPA-capable molecule, and
$c$ is some function of system parameters. The function $p$ may be zero-order
in its arguments, in which case RPA is ‘unconstrained’ by any other system
properties. In all other cases, Equation 2 encodes a concept known as
constrained integral control.

In simple biochemical network designs, a linear change of coordinates may be
sufficient to identify the internal model, such that an integrator of the form
in Equation 2 lies in the rowspan of the system. Antithetic integral control
[16] represents a well-known example of this type and is, topologically
speaking, a feedback-regulated system. A simple representation of this control
motif is depicted in Figure 1A, with four interacting molecules with
concentrations $A,B,C$ and $D$. From the reaction forms shown, is it clear
that

$\frac{dC}{dt}-\frac{dD}{dt}=k_{4}\left(B-\frac{k_{6}}{k_{4}}\right),$

which captures the internal model explicitly, since $B^{*}=k_{6}/k_{4}$
represents the setpoint of the molecule $B$. New methods have recently been
developed to detect an internal model in the rowspan of an RPA (or
ACR)-capable system, if such exists, in more complex networks of biochemical
interactions [22, 27].

Figure 1: Two simple motifs from which an internal model may be identified
through a simple change of coordinates. (A) An antithetic integral control
motif [16], as discussed by [28]. (B) An incoherent feedforward motif [26].

Feedforward structured networks can also be recast as integral feedback
systems through a judicious (generally non-linear) change of coordinates. Such
a nonlinear mapping may be difficult to identify for complex feedforward-
structured networks, but may be readily obtained for particularly simple
feedforward motifs. For instance, Shoval et al. [29] consider a minimal
incoherent feedforward motif of the form

$\displaystyle\frac{dx_{1}}{dt}$ $\displaystyle=\alpha u-\delta x_{1},$
$\displaystyle\frac{dx_{2}}{dt}$ $\displaystyle=\beta\frac{u}{x_{1}}-\gamma
x_{2},$

as depicted in Figure 1(B), and show that the map

$\left(X_{1},X_{2}\right)\to\left(z_{1},z_{2}\right)=(X_{2},\alpha
X_{2}-\beta\textrm{log}X_{1})$

is a diffeomorphism with inverse $X_{1}=e^{(\alpha z_{1}-z_{2})/\beta}$ and
$X_{2}=z_{1}$, which transforms the reaction equations into the partitioned
form

$\displaystyle\frac{dz_{1}}{dt}$ $\displaystyle=\beta ue^{(z_{2}-\alpha
z_{1})/\beta}-\gamma z_{1},$ $\displaystyle\frac{dz_{2}}{dt}$
$\displaystyle=\beta\delta-\alpha\gamma z_{1}.$

The second of these equations captures the internal model, since the solutions
to $dz_{2}/dt=0$ represent all possible constant signals, and determine the
setpoint of the system to be $z_{1}=X_{2}=\delta\beta/{\alpha\gamma}$.

The recognition that RPA may be achieved even for very complex networks by
adding an integrator to the network within a feedback structure has been of
tremendous utility for synthetic implementations [8]. On the other hand, this
recognition alone does not provide a straightforward tool for exploring a more
comprehensive design space of RPA-capable networks - particularly in view of
the fact that the connection to integral feedback regulation may be very
complex for networks without a clear (physical) feedback structure (such as
feedforward structured networks) [26]. Moreover, there is now strong
experimental evidence that RPA in natural (endogenous) cellular chemical
reaction networks frequently requires multiple adaptation mechanisms working
together in concert [30]. Even in the relatively simple chemotaxis-regulating
cascades in organisms such as Dictyostelium, for example, it appears that
multiple distinct RPA modules are at work [30], making it difficult to propose
synthetic implementations of these kinds of network structures relying only on
the concept of integral control.

In the next section, we consider a complementary approach to the integral
control viewpoint, by considering how an algebraic representation of the RPA
problem may be analysed in such a way as to reveal the complete, comprehensive
design space of RPA-capable network topologies. These methods identify RPA-
capable network structures directly, and retain a representation of the RPA
problem in terms of the original coordinates (ie. the concentrations of the
interacting molecules). Rather than extract a single internal model for the
constant input stimuli to an RPA-capable system, this approach reveals how
such an internal model could be distributed topologically over potentially
vast and intricate networks of interacting molecules.

###  III-B The topological RPA problem

Many algebraic formulations for the RPA problem have been proposed [3, 31,
32], and have been shown to be equivalent to the requirement for integral
control [3, 31]. Commonly, these formulations allow the network to have
multiple ‘input nodes’, ie. molecules that are subjected to stimulation or
perturbation by an external source, as well as multiple ‘output nodes’, which
are the molecular endpoints of interest, which are expected to exhibit the RPA
property. But to delineate the full solution space of RPA-capable network
topologies from such an algebraic condition, it is expedient to start from a
simpler assumption: a network of interacting nodes that are perturbed at a
single input node, and for which a single output node is monitored for the RPA
property, thereby allowing network topologies to be induced relative to the
designated input-output node pair [3]. For an RPA-capable network topology
identified by this method, there may exist additional nodes to which a
stimulus may be applied while preserving the RPA property at the designated
output node. Likewise, an RPA-capable network design may also exhibit RPA at
nodes other than the designated output node. In this way, a complete picture
of RPA-capable multi-input/multi-output network designs emerges.

The topological RPA problem has now been solved in complete generality [3],
and is the fulfilment of a longstanding conjecture in biological network
theory [33, 34, 35, 36, 37] which held that robust biosystems are likely to
require very specific network architectures - a conjecture based largely on
analogies with engineering control systems and the well-defined design
features that are known to be required, or at least beneficial, in those
contexts. The essential ideas underlying this general solution are as follows:

Suppose we have a network consisting of $n$ interacting molecules, or nodes,
denoted $N_{1},N_{2},\ldots,N_{n}$. One of these nodes ($N_{i}$, say) will be
selected as the designated input node, and one node ($N_{j}$, not necessarily
distinct from $N_{i}$) will be selected as the output. To each non-input node,
we associate a reaction rate $dN_{k}/dt=f_{k}(N_{1},\ldots,N_{n})$, which
encodes both the strength and the nature (stimulating or inhibitory) of the
inter-molecular interactions upon the node $N_{k}$. The input node is also
regulated by an external stimulus, $S$, and has a reaction rate
$dN_{i}/dt=f_{i}(N_{1},\ldots,N_{n},S)$. From this simple mathematical
representation of the interactions in an $n-$node network, and its encounter
with an external stimulus, it is straightforward to show that the network is
capable of RPA at node $N_{j}$ only if

$\frac{\textrm{Det}(M_{ij})}{\textrm{Det}(J_{n})}=0,$ (3)

where
$J_{n}=\frac{\partial(f_{1},\ldots,f_{n})}{\partial(N_{1},\ldots,N_{n})}$
denotes the system Jacobian, and $M_{ij}$ is the Jacobian minor obtained by
eliminating the input row and the output column from $J_{n}$. Thus, a
necessary condition for RPA in any network, regardless of size or complexity,
is

$\textrm{Det}(M_{ij})=0,$ (4)

with

$\textrm{Det}(J_{n})\neq 0.$ (5)

Equation 4 is known as the RPA equation [3], which must be satisfied in any
RPA network, subject to the RPA constraint (Equation 5).

We note that the RPA equation will, in general, contain a very large number of
terms, depending on the number of nodes and the interconnectedness of the
network, highlighting the fact that the topological RPA problem quickly eludes
computational screening approaches as $n$ grows. Indeed, for a fully-connected
network, where every node regulates every other node, there will be $(n-1)!$
terms in the RPA equation, each containing $(n-1)$ factors. We illustrate in
Figure 2 how quickly the size of the RPA equation grows with the number of
interacting nodes in the network.

Figure 2: The RPA Equation (Equation 4) grows factorially with the number of
interating nodes, $n$ in the network.

The cornerstone of the methodology developed in [3] is the interpretation of
the determinant of $M_{ij}$ not simply as a function Det:
$\mathbb{R}^{(n-1)\textrm{x}(n-1)}\to\mathbb{R}$, but as a set of $\tau$
terms, $1\leq\tau\leq(n-1)!$, with each term carrying either a positive or a
negative sign. Crucially, the $(n-1)$ factors of each term correspond in a
very precise way to topological characteristics of the underlying network,
representing either structural features of the network, or constraints on
reaction kinetics (see Figure 3).

Figure 3: Relationship of the factors within individual terms of the RPA
Equation to the topology of the underlying network.

It is clear from this set theoretic interpretation of the RPA Equation that
partitions of the set may exist for which the contents of every subset can sum
to zero independently of the contents of every other subset [3]. By
identifying general conditions under which such ‘independently adapting
subsets’ can exist, we can determine, definitively, the topological structure
that any RPA-capable network must possess, no matter how large or complex the
network. In particular, it has been shown [3] that there exist two and only
two possible types of ‘independently adapting subset’ of the RPA equation,
known as S-sets and M-sets. S-sets contain terms that are all identically
zero, for all possible values of input stimulus, and all possible parameter
values, and arise when a feedback structure exists within the associated
network, and where one or more nodes embedded into that feedback structure
satisfy a special constraint on reaction kinetics known as Opposer kinetics.
These S-sets correspond to the presence of a special subnetwork structure
known as an Opposer Module (see Figure 4(A)). By contrast, M-sets contain
terms that are all strictly non-zero, but together sum to zero for all
possible values of input stimulus, and all possible parameter values. M-sets
arise when a collection of parallel pathways exists in the network, and where
the nodes embedded into those parallel pathways satisfy a special constraint
on reaction kinetics known as Balancer kinetics. In addition, the convergence
point for these parallel pathways must also satisfy a constraint on reaction
kinetics known as Connector kinetics in order to complete the ‘balancing act’
of the M-set. M-sets correspond to the presence of a special subnetwork
structure known as a Balancer Module (see Figure 4(B)).

Figure 4: Two module classes represent a topological basis for all RPA-capable
networks. RPA-capable nodes are indicated by a red asterisk; all other nodes
are unable to exhibit RPA. $S_{N}$ indicates that an arbitrary subnetwork may
be embedded at the location indicated without affecting the RPA capacity of
the network’s output node. A complete description of RPA basis modules is
given in [3]. (A) An Opposer Module, which involves a feedback structure. At
least one opposer node, which satisfies a special constraint known as opposer
kinetics, must be embedded into the feedback portion of an Opposer Module. In
the representation depicted here, two opposer nodes (indicated in yellow)
collaborate to form a special architecture known as a two-node opposing set.
(B) A Balancer Module is characterised by a collection of parallel pathways,
with a diverter node (D) at the source point of the parallel pathways, and a
connector node (C) at the downstream convergence point. All nodes between D
and C, known as balancer nodes (indicated in blue), must satisfy a special
constraint known as balancer kinetics. The connector node must also satisfy a
special constraint known as connector kinetics.

Opposer and Balancer Modules are generalizations of the simple three-node RPA-
capable motifs identified by Ma et al. [38] through computational screening
methods. But these more general RPA-promoting modules may be vastly more
complex, and can contain previously unknown topological features such as
Opposing sets. We refer interested readers to [3] and [39] for a comprehensive
description of Opposer and Balancer Modules.

These two RPA Modules thus represent a topological basis for all RPA-capable
networks. In other words, RPA-capable networks are necessarily modular, and
must be decomposable into these well-defined classes of network substructures
- an intriguing feature of immense practical utility in the context of
synthetic network design.

Diagrammatic methods have recently been developed to assist in the analysis of
RPA capacity for specific network architectures [39]. For the network
structure given in Figure 5, for instance, it can easily be shown that the RPA
equation for this particular network contains just three terms, and admits two
(and only two) decompositions into RPA basis modules: a single Opposer Module
(Figure 5(A)), or a Balancer Module in parallel with a (smaller) Opposer
Module (Figure 5(B)). The rules of interconnectivity for RPA basis modules to
create larger, multi-modular networks, and the concept of series and parallel
connections of basis modules, are also known in complete generality, and are
explained in full in [3] and [39].

Figure 5: A simple network architecture, which supports two different
decompositions into RPA basis modules. Opposer nodes are indicated in yellow;
balancer nodes are indicated in blue; and connector nodes are indicated in
green. It can readily be shown via diagrammatic methods [39] that the RPA
Equation for this network contains three terms, of which two are positive and
one is negative. These terms admit two different partitions into
‘independently adapting subsets’, as indicated schematically beneath each
network. (A) An Opposer Module emerges as the three terms of the RPA equation
are assigned to a single S-set (yellow). (B) An Opposer Module in series with
a Balancer Module emerges as the negative term and one of the positive terms
are assigned to an M-set (blue), while the remaining term is assigned to an
S-set (yellow).

As we close this section, we note that the topological view of the RPA problem
brings into sharp relief the concept that biochemical signalling networks do
not simply transmit biochemical signals, but are required to perform precise
and coordinated computations on those signals in order to implement robust
function. Indeed, unlike the well-established control theory perspective,
which emphasizes the requirement for an internal model (and integral feedback
in the case of RPA), topological analysis of robust network designs highlights
the presence of a constellation of invariants, which includes the key
invariant referred to as the system setpoint, which are distributed in very
specific topological arrangements across a potentially large and intricate
signalling network. Figure 6 illustrates this key idea via a schematic
representation of invariant distributions in two minimal RPA modules - a
three-node Balancer Module (Figure 6(A)), and an Opposer Module containing a
topological feature known as a three-node opposing set (Figure 6(B)). Each
node that is subjected to a constraint on its reaction kinetics - a balancer
node, an opposer node or a connector node - thereby computes an RPA-promoting
invariant in order to confer the RPA property on the output node of the
module, and ultimately the entire network.

Figure 6: RPA network topologies compute a collection of invariants which,
together, confer the RPA property on the output node. (A) A minimal Balancer
Module, in which balancer kinetics at node B compute a balancer invariant,
which is then transferred to node C, which then computes the setpoint. (B) A
three-node opposing set. Each of the three independent opposer nodes (yellow)
calculates its own opposer invariant, which is transferred from $O_{1}$ to
$O_{2}$ to $O_{3}$, thereby conferring RPA on $P_{2}$.

##  IV Bayesian Perspectives on Design

Bayesian inference procedures [40] arguably provides a very natural framework
for expressing and reporting uncertainty. Because of this, and its ability to
accomodate background information and knowledge it has become an increasingly
popular in systems biology [41]. One of the main advantages is that it
provides a powerful and internally consistent framework for model selection –
that is, choosing from a set of candidate models, those which have a high
probability or producing the observed data [42]. This can also be used in the
context of design for synthetic biology.

###  IV-A Bayesian Inference for Design

Bayesian design approaches take a very different, but we feel, complementary
approach. Here we consider sets of models
$\mathcal{F}=\\{f_{1},f_{2},\ldots,f_{N}\\}$, and try to identify which of
these models is best capable of meeting the design objectives, $T$. Here we
adopt an inferential perspective, but instead of calibrating a model against
observed data we calibrate the models against the data we would like to see,
that is, the design objective.

For each model $f_{i}\in\mathcal{F}$ the aim is to obtain an estimate for the
probability that $f_{i}$ is the best model capable of producing the target
behaviour, given $T$ and the set of candidate models, $\mathcal{F}$,
$\text{Pr}(f_{i}|T)$. This posterior model probability is given by

$\text{Pr}(f_{i}|T)=\frac{\text{Pr}(T|f_{i})\pi(f_{i})}{\text{Pr}(T)}=\frac{\text{Pr}(T|f_{i})\pi(f_{i})}{\sum_{i}\text{Pr}(T|f_{i})\pi(f_{i})},$
(6)

where $\pi(f_{i})$ is the prior probability of model $f_{i}$,
$\text{Pr}(f_{i})$, denotes the likelihood, and
$\text{Pr}(T)=\sum_{i}\text{Pr}(T|f_{i})\pi(f_{i})$, the evidence, which can
be interpreted as a normalising factor (akin to the partition function in
statistical physics).

The terms in Equation (6) refer to the models. We have subsumed the parameter
dependence implicitly into these terms; but we need to consider them, of
course, as a model $f_{i}$ will typically only exhibit the target behaviour,
$T$, for certain regions of the parameter spaces, $\Omega_{\theta_{i}}$. We
let the $i$ subscript refer to the model, $f_{i}$, and we allow different
parameter spaces for the different models. The marginal likelihood,
$\text{Pr}(T|f_{i})$, is given by

$\text{Pr}(T|f_{i})=\int_{\Omega_{i}}d\theta^{\prime}\text{Pr}(T|f_{i},\theta^{\prime})\pi_{i}(\theta^{\prime}),$
(7)

where $\text{Pr}(T|f_{i},\theta)$ is the likelihood of the parameter, and
$\pi_{i}(\theta)$ denotes the prior of the parameters for model $i$
(accounting for the fact that parameter vectors can differ between models).

The rationale for Bayesian design is that it allows us to estimate the
probability that a given model, $f_{i}$ will be able to meet the design
objectives, $T$. This is dependent not only on $T$, but also on the set of
alternative models considered. The model posterior probability,
$\text{Pr}(f_{i}|T)$, can thus be used to choose from the available
alternative designs.

This approach has an intrinsic robustness criterion embedded at its heart: in
Equation (7) the integral is taken over the parameter space, which means that
the model, $f_{i}$, with the maximum a posteriori probability will generally
show good performance against the design objectives over a range of
parameters. With this information it is possible to assess the likely
robustness – a concept that we will revisit below – of a design, and compare
competing designs in terms of expected performance and robustness.

###  IV-B Approximate Bayesian Computation for Design

This likelihood is also difficult to evaluate for most important problems of
scientific and technological interest. As a result, a range of approximations
have been developed that forgo the explicit evaluation of the likelihood in
favour of using simulations. Of these, approximate Bayesian computation [43,
44], or ABC, is perhaps the best known alternative to evaluating the full
posterior via the likelihood function. It has been used in the context of
design with some success[4] and we outline it below before considering
applications, and generalisations.

In brief, we sample a parameter value, $\theta^{\prime}$ from the prior,
$\pi_{i}(\theta)$, and we then simulate the model for that parameter value to
obtain simulated data $\mathcal{D}=\\{d_{1},d_{2},\ldots,d_{n}\\}$. This is
then compared to the real data – or in the present context to the design
objectives, perhaps via surrogate data $\mathcal{D}_{T}$ – via a distance
function, $\Delta(X,Y)$, such that $\theta^{\prime}$ is regarded as a draw
from the approximate posterior, if

$\Delta(\mathcal{D},T)\leq\epsilon$ (8)

where $\epsilon$ has to be sufficiently small. The full posterior for the
parameter of a given model is thus replaced by

$\text{Pr}(\theta|T)\approx\text{Pr}_{\text{ABC}}(\theta|\Delta(\mathcal{D},T)\leq\epsilon),$
(9)

and analogously for the model posterior. (Here we can think of sampling models
and then model parameters simultaneously). Note that here we do not, as is
frequently the case[45], take a summary statistic of real and simulated data,
but compare the simulation output directly to the data (or the design
objectives).

In ABC inference the choice of the distance function is of pivotal importance
and can make profound differences on the quality and indeed validity of the
inference. Classical rejection-based ABC (where only a single threshold is
chosen) tends to be computationally inefficient and wasteful. Instead most
inference approaches in the literature adopt a sequential Monte Carlo (SMC)
approach, where a sequence of thresholds,
$\epsilon_{1}\geq\epsilon_{2}\geq\ldots\geq\epsilon_{n}$ is used (either
predefined or constructed “online” over the course of the inference procedure
– thus either $n$ or $\epsilon_{n}$ can be a random variable). Choice of this
threshold schedule and other aspects of the ABC-SMC approach are still
requiring a lot of attention.

In the context of design the distance function also needs to do a lot of heavy
lifting, both in terms of performance and, more fundamentally, for the
conceptual set-up of the ABC design approach.

In the simplest case we could take an indicator function,

$\Delta(\mathcal{D},T)=\begin{cases}0\qquad\text{ iff $\mathcal{D}$ has
property $T$}\\\ 0\qquad\text{ otherwise.}\end{cases}$ (10)

But if we want to put ABC design to work in the context of the vastly more
efficient SMC samplers we have to rethink such indicators. This may not always
be possible, of course. The apparent simplicity of ABC approaches masks that
there is a great deal of care required to make sure that model selection can
be performed in this framework (discussed e.g. in [45, 46, 47]

###  IV-C ABC Design for Perfect Adaptive Behaviour

We can illustrate the ABC for design approach in the context of adaptive
behaviour. In addition to the algebraic methods [3], which can be used to
distill the design characteristics, large-scale computational searches [38]
have been used. These identified two models that are able to generate adaptive
behaviour. Either the algebraic approach or brute-force computational surveys
can be used to produce sets of model designs,
$\mathcal{F}=\\{f_{1},\ldots,f_{n}\\}$ which can then be used as the starting
point for model selection for design.

The output characteristics in [4] were chosen to be the adaptation efficiency
and sensitivity [38], given by

$\displaystyle E=$
$\displaystyle\frac{(O_{2}-O_{1})/O_{1}}{(I_{2}-I_{1})/I_{1}}$ (11) and
$\displaystyle S$
$\displaystyle\frac{(O_{\text{peak}}-O_{1})/O_{1}}{(I_{2}-I_{1})/I_{1}},$ (12)

where $I_{1}$ and $I_{2}$ are the baseline and activated input levels; $O_{2}$
and $O_{1}$ are the steady state outputs before and after the signal is
applied; and $O_{\text{peak}}$ is the maximal transient output level.

In Fig. LABEL:fig:fig_abc we show two models (A), (B), that have been argued
to exhibit adaptive behaviour in previous studies [4, 38], and their
respective posterior probabilities to produce adaptive behaviour if
cooperativity is allowed (C) or not allowed (D). The two posterior
probabilities are sufficiently different in the cooperative scenario (C) to
warrant choice of model 1 $\text{Pr}(f_{1}|T)\approx 0.96$; whereas for non-
cooperative dynamics the choice is less clear with $\text{Pr}(f_{2}|T)\approx
0.64$. This is not decisive enough and further analysis, and consideration of
further design candidates, may become necessary.

##  V Robustness in Synthetic Biology Design

Below we develop two notions of robustness in synthetic biology design:
Biological Robustness refers to the ability of a system to perform its
function in the presence of perturbations.
Model Robustness is related to mismatches between the mathematical model and
the real biological systems it sets out to describe.
Both are relevant in the context of synthetic and engineering biology and
design approaches, and are discussed below.

###  V-A Robustness of Biological Systems and Models

Robustness as a property of both biological systems and their representations
as mathematical models has been the subject of much discussion [11, 12, 48,
49, 50, 51, 52, 53, 54, 55, 56, 57]. For biological systems, robustness is
generally defined as the property that ensures the system maintains its
functions in the presence of perturbations from internal and external sources
[33, 58, 34]. In the context of design we take a pragmatic point of view of
robustness. We are interested in identifying designs that are robust to (i)
inherent stochasticity; (ii) environmental and physiological signals and
fluctuations; (iii) uncertainty in our knowledge of the precise system
architecture and kinetic reaction rates.

For mathematical models many working definitions of robustness have been
proposed, which are based on single- and multi-parameter sensitivity analysis
[50, 57]. They include, for example, invariance under external fluctuations
affecting the initial conditions [48]; or structural stability [59]. There is,
however, no clear and consistent definition that allows for universal and
consistent quantification of robustness, which is necessary if the robustness
characteristics of different systems are to be compared. Perhaps the primary
importance of quantifying model robustness is in model validation, as a valid
model must exhibit equivalent robustness to the corresponding biological
system it represents. In fact, a lack of robustness may indicate a weak model
that provides a poor representation of the biological system [51], so
robustness may be used as a selection tool for promising models that correctly
represent the functions of the biological system. The determination of weak
models through robustness quantification is also useful for providing insight
into the robustness of the biological system in a synthetic biology design
context.

Importantly, robustness relates to the maintenance of the functions of the
system and not necessarily the states of the system, so is more general than
stability and homeostasis and may involve transitions to new steady states or
even to an instability if required to maintain function [34]. For example,
system function where a morphogen concentration gradient determines the
position in a biological system appears to be robust despite temporal
perturbations of the gradient to different steady-states [60]. Despite their
complexity, evolved biological systems apparantly exhibit robustness as a
ubiquitous and fundamental property [33, 61, 62]. Synthetic cellular systems
have to work in the context of such naturally evolved contexts.

###  V-B Mathematical Definitions of Robustness

Unfortunately, and as noted some time ago [34], there is still no adequate
mathematical theory of robustness. There are a number of analytical issues to
be addressed, in particular the difficulty involved with high-dimensional
parameter spaces. For example, control-theoretic techniques have been employed
to study the robustness of biochemical network models [50]; however, these
approaches are limited to low-dimensional parameter spaces [33].

A robust model must maintain the represented biological function with respect
to some degree of perturbation of the inputs. In other words, we require small
changes in the inputs – model parameters, intitial conditions, boundary
conditions, external conditions – to produce sufficiently small changes in the
output so that the integrity of the biological function remains. Note that
this is different from requiring invariance of the output: the system output
can vary as long as the biological function is maintained. We therefore only
require invariance of biological function: the design objectives, $T$ still
have to be met.

We require a way of assessing robustness and focus here, for didactic reasons,
on robustness to parametric uncertainty or variation. In order to compare
distances between input values and output values (including any design
objective, $T$) we require the sets of input and output values to form metric
spaces.

We furthermore require a functional relationship between the input and output
metric spaces, which would generally be a continuous function in order to
ensure that the output values change smoothly with changes in the inputs. It
is also possible that an appropriate function has discontinuities, such as a
step discontinuity, so may be only piecewise continuous. In this case we just
need to ensure that the function output still corresponds to invariance of
biological function by accounting for the behaviour near the discontinuities.
This makes it possible to accommodate qualitative design criteria. We would
also like to be able to specify the robustness of a system by a single
parameter; we therefore need a methodology that can be applied to the whole
parameter space to give an overall measure of robustness.

From a global perspective, we require that close values of the input produce
output values that are consistently close. Note that local continuity of a
function, or global uniform continuity, are not sufficient for defining
robustness as they allow for arbitrarily large derivatives of the function,
and thus a lack of robustness. The Lipschitz function [63, Section 9.4, page
154] immediately lends itself to a rigorous definition of robustness for
mathematical models. But inexplicably it has not been employed within the
mathematics literature, to our knowledge; it has, however, appeared in other
fields such as software engineering [64] and cyber-physical systems [65].

We denote the input and outputs by $X$ and $Y$, respectively; and we define
suitable metrics, $d_{X}$ and $d_{Y}$ on the statespaces, for $X$ and $Y$; we
can therefore consider $(X,d_{X})$ and $(Y,d_{Y})$ to be metric spaces. A
function $f\colon X\to Y$ is a _Lipschitz function_ on $X$ with a _Lipschitz
constant_ $K\geq 0$ if $f$ satisfies

$d_{Y}(f(x_{1}),f(x_{2}))\leq K\,d_{X}(x_{1},x_{2})$ (13)

for all $x_{1}$, $x_{2}\in X$. We also say that $f$ is _Lipschitz continuous_.
The smallest such constant $K$ is _the_ Lipschitz constant. The Lipschitz
constant provides an upper bound on the change of the function output given a
change in the input. Ideally we would like an analytical value for the
Lipschitz constant, however there may be Lipschitz functions for which this is
not easily obtained, so numerical estimates may be helpful. Numerical methods
for estimating Lipschitz constants exist [66]; however, there is no guarantee
that the estimated constant is valid globally. In control theory _transfer
functions_ are employed to provide a functional relationship between the
inputs and outputs of a system, often within the frequency domain, by applying
the Laplace transform to both the inputs and outputs, thereby allowing direct
comparison of functions [67]. While the overall approach using transfer
functions overlaps with our requirements, it is perhaps too specific and
limited for characterising the robustness of complex biological systems.

###  V-C Parametric Robustness as a Design Criterium in Synthetic Biology

The algebraic approaches outlined above allow us to detect the design
principles underlying a given desired behaviour or objective, $T$. This allows
us to design models that meet the objectives by construction.

The Bayesian approach contains an assessment of robustness via the
posterior[68], which integrates over the parameters, Equation 7. Rather than
the design objective, $T$, we may use the posterior parameter density,
$\text{Pr}(\theta|T)$, or the ABC posterior, (7), as the output, and then
compare $d_{\theta}(\theta,\theta^{\prime})$ and
$d_{P}(\text{Pr}(\theta|T),\text{Pr}(\theta^{\prime}|T)$ using the Lipschitz
criterion (13). We can focus this analysis on, for example, the 90% credible
interval [40].

If our design, that is the mathematical model, is robust against variation in
parameters, if it is structurally stable, and if is is even robust against
changes to the model structure [69], then we can view this as reassuring that
the design will also be robust to the discrepancies between mathematical model
and biological, reality[41, 6].

##  VI Discussion and Conclusions

The models in section III and in Figure 6 are structurally very similar to the
models considered in section V and Figure LABEL:fig:fig_abc, but differ in
practical terms as we first use mass action kinetics and then Michaelis-Menten
kinetics. Algebraic methods are more powerful wherever they are applicable, as
they result in definitive criteria that must be satisfied for a given
behaviour to emerge. From this we can often gain profound insights into the
dynamics of cellular systems [3, 5].

For many classes of design objective, other methods need to be applied, and
the Bayesian design framework is one method that is widely applicable for
different models. Basing design on Bayesian model selection [70, 71, 4] has
two intrinsic advantages: (i) we are able to balance the performance and the
robustness of different designs; and (ii) for ABC methods we can, at least in
principle, apply Bayesian design to arbitrary models as long as we are able to
simulate them. The first point is potentially an advantage over large-scale or
exhaustive analyses of models [38, 72, 69, 6]. There, robustness can be
assessed post-hoc. While global analyses and Bayesian methods are both
computationally expensive, the former are better suited for the comprehensive
analysis of a large number of systems: the objective there is to scan model
and parameter spaces to find if there are any parameters for which the design
objectives, $T$, are realised. In the Bayesian approach, by contrast, we also
obtain estimates of the parameters, $\theta_{i}$, for each model
$f_{i}\in\mathcal{F}$.

So we can view the algebraic and Bayesian frameworks as complementary, and
given the complexities of rationally designing synthetic circuits [1, 73] a
flexible and diverse tool-box of alternative design methods is clearly
beneficial. But algebraic methods can also support and work in conjunction
with Bayesian methods (or exhaustive scans of large sets of models): we can
determine initial sets of models, $\mathcal{F}=\\{f_{1},\ldots,f_{n}\\}$ to
fulfil the design criteria based on algebraic rules or design principles (even
if under perhaps restricted kinetic regimes, such as mass-action kinetics); a
natural choice for the prior model probabilities would then be,
$\pi(f_{i})=1/n$ for all $i=1,\ldots,n$. These models can then be subjected to
the Bayesian model selection/design process to choose those models that have
the highest probability of delivering the design objective, $T$ (and do so
robustly).

The challenge of design in synthetic biology is considerable. Unlike in
physics and many engineering disciplines we cannot draw on first-principles.
Reliable mathematical models for cellular systems are still, despite much
recent progress, rare, but we need them to apply rational engineering design
approaches in synthetic biology. The methods described here – algebraic ways
to identify design principles; Bayesian design; and the concept of robustness
– can help to identify synthetic circuits capable of delivering a target
objective, $T$. In combination they allow us to take some of the guess-work
and trial-and-error out of synthetic biology design.

## References

  * [1] W. A. Lim, “Designing customized cell signalling circuits.” _Nature Reviews Molecular Cell Biology_ , vol. 11, no. 6, pp. 393–403, 2010.
  * [2] J. Santos-Moreno, E. Tasiudi, J. o. r. Stelling, and Y. Schaerli, “Multistable and dynamic crispri-based synthetic circuits.” _Nature communications_ , vol. 11, no. 1, pp. 2746–2748, 2020.
  * [3] R. P. Araujo and L. A. Liotta, “The topological requirements for robust perfect adaptation in networks of any size.” _Nature Communications_ , vol. 9, no. 1, p. 1757, 2018.
  * [4] C. P. Barnes, D. Silk, X. Sheng, M. P. H. Stumpf, and M. P. H. Stumpf, “Bayesian design of synthetic biological systems.” _Proceedings of the National Academy of Sciences of the United States of America_ , vol. 108, no. 37, pp. 15 190–15 195, 2011.
  * [5] H. A. Harrington, E. Feliu, C. Wiuf, M. P. H. Stumpf, and M. P. H. Stumpf, “Cellular compartments cause multistability and allow cells to process more information.” _Biophysical journal_ , vol. 104, no. 8, pp. 1824–1831, 2013\.
  * [6] N. S. Scholes, D. Schnoerr, M. Isalan, and M. P. H. Stumpf, “A comprehensive network atlas reveals that turing patterns are common but not robust.” _Cell Systems_ , vol. 9, no. 3, pp. 243–257.e4, 2019.
  * [7] D. Muzzey, C. A. Gómez-Uribe, J. T. Mettetal, and A. van Oudenaarden, “A systems-level analysis of perfect adaptation in yeast osmoregulation,” _Cell_ , vol. 138, no. 1, pp. 160–171, 2009.
  * [8] M. H. Khammash, “Perfect adaptation in biology,” _Cell Systems_ , vol. 12, no. 6, pp. 509–521, 2021.
  * [9] U. B. Kaupp, “Olfactory signalling in vertebrates and insects: differences and commonalities,” _Nature Reviews Neuroscience_ , vol. 11, no. 3, pp. 188–200, 2010.
  * [10] K.-W. Yau and R. C. Hardie, “Phototransduction motifs and variations,” _Cell_ , vol. 139, no. 2, pp. 246–264, 2009.
  * [11] N. Barkai and S. Leibler, “Robustness in simple biochemical networks,” _Nature_ , vol. 387, pp. 913–917, 1997.
  * [12] U. Alon, M. G. Surette, N. Barkai, and S. Leibler, “Robustness in bacterial chemotaxis,” _Nature_ , vol. 397, pp. 168–171, 1999.
  * [13] T. Toni, G. Jovanovic, M. Huvet, M. Buck, and M. P. H. Stumpf, “From qualitative data to quantitative models: analysis of the phage shock protein stress response in escherichia coli.” _Bmc Systems Biology_ , vol. 5, no. 1, p. 69, 2011.
  * [14] N. Joly, C. Engl, G. Jovanovic, G. Jovanovic, M. Huvet, T. Toni, X. Sheng, M. P. H. Stumpf, M. P. H. Stumpf, and M. Buck, “Managing membrane stress: the phage shock protein (psp) response, from molecular mechanisms to physiology.” _FEMS microbiology reviews_ , vol. 34, no. 5, pp. 797–827, 2010\.
  * [15] T. You, P. Ingram, M. D. Jacobsen, E. Cook, A. McDonagh, T. Thorne, M. D. Lenardon, A. P. De Moura, M. C. Romano, M. Thiel _et al._ , “A systems biology analysis of long and short-term memories of osmotic stress adaptation in fungi,” _BMC research notes_ , vol. 5, no. 1, pp. 1–16, 2012.
  * [16] C. Briat, A. Gupta, and M. Khammash, “Antithetic integral feedback ensures robust perfect adaptation in noisy biomolecular networks,” _Cell systems_ , vol. 2, no. 1, pp. 15–26, 2016.
  * [17] F. Mairet, “A biomolecular proportional integral controller based on feedback regulations of protein level and activity,” _Royal Society open science_ , vol. 5, no. 2, p. 171966, 2018.
  * [18] X. Y. Ni, T. Drengstig, and P. Ruoff, “The control of the controller: molecular mechanisms for robust perfect adaptation and temperature compensation,” _Biophysical journal_ , vol. 97, no. 5, pp. 1244–1253, 2009\.
  * [19] A. Eldar, R. Dorfman, D. Weiss, H. Ashe, B.-Z. Shilo, and N. Barkai, “Robustness of the bmp morphogen gradient in drosophila embryonic patterning,” _Nature_ , vol. 419, no. 6904, pp. 304–308, 2002.
  * [20] D. Ben-Zvi and N. Barkai, “Scaling of morphogen gradients by an expansion-repression integral feedback control,” _Proceedings of the National Academy of Sciences_ , vol. 107, no. 15, pp. 6924–6929, 2010.
  * [21] F. Xiao and J. C. Doyle, “Robust perfect adaptation in biomolecular reaction networks,” in _2018 IEEE Conference on Decision and Control (CDC)_. IEEE, 2018, pp. 4345–4352.
  * [22] D. Cappelletti, A. Gupta, and M. Khammash, “A hidden integral structure endows absolute concentration robust systems with resilience to dynamical concentration disturbances,” _Journal of the Royal Society Interface_ , vol. 17, no. 171, p. 20200437, 2020.
  * [23] C. I. Byrnes and A. Isidori, “Limit sets, zero dynamics, and internal models in the problem of nonlinear output regulation,” _IEEE Transactions on Automatic Control_ , vol. 48, no. 10, pp. 1712–1723, 2003.
  * [24] B. A. Francis and W. M. Wonham, “The internal model principle for linear multivariable regulators,” _Applied mathematics and optimization_ , vol. 2, no. 2, pp. 170–194, 1975.
  * [25] ——, “The internal model principle of control theory,” _Automatica_ , vol. 12, no. 5, pp. 457–465, 1976.
  * [26] E. D. Sontag, “Adaptation and regulation with signal detection implies internal model,” _Systems & control letters_, vol. 50, no. 2, pp. 119–126, 2003.
  * [27] R. L. Karp, M. P. Millán, T. Dasgupta, A. Dickenstein, and J. Gunawardena, “Complex-linear invariants of biochemical networks,” _Journal of theoretical biology_ , vol. 311, pp. 130–138, 2012.
  * [28] J. E. Ferrell Jr, “Perfect and near-perfect adaptation in cell signaling,” _Cell systems_ , vol. 2, no. 2, pp. 62–67, 2016.
  * [29] O. Shoval, U. Alon, and E. Sontag, “Symmetry invariance for adapting biological systems,” _SIAM Journal on Applied Dynamical Systems_ , vol. 10, no. 3, pp. 857–886, 2011.
  * [30] O. Hoeller, D. Gong, and O. D. Weiner, “How to understand and outwit adaptation,” _Developmental cell_ , vol. 28, no. 6, pp. 607–616, 2014.
  * [31] T.-M. Yi, Y. Huang, M. I. Simon, and J. Doyle, “Robust perfect adaptation in bacterial chemotaxis through integral feedback control,” _Proceedings of the National Academy of Sciences_ , vol. 97, no. 9, pp. 4649–4653, 2000.
  * [32] Z. F. Tang and D. R. McMillen, “Design principles for the analysis and construction of robustly homeostatic biological networks,” _Journal of theoretical biology_ , vol. 408, pp. 274–289, 2016.
  * [33] H. Kitano, “Biological robustness,” _Nature Reviews Genetics_ , vol. 5, pp. 826–837, 2004.
  * [34] ——, “Towards a theory of biological robustness,” _Molecular Systems Biology_ , vol. 3, p. 137, 2007.
  * [35] R. P. Araujo, L. A. Liotta, and E. F. Petricoin, “Proteins, drug targets and the mechanisms they control: the simple truth about complex networks,” _Nature Reviews Drug discovery_ , vol. 6, no. 11, pp. 871–880, 2007.
  * [36] R. P. Araujo and L. A. Liotta, “A control theoretic paradigm for cell signaling networks: a simple complexity for a sensitive robustness,” _Current opinion in chemical biology_ , vol. 10, no. 1, pp. 81–87, 2006.
  * [37] J. J. Tyson, K. C. Chen, and B. Novak, “Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell,” _Current opinion in cell biology_ , vol. 15, no. 2, pp. 221–231, 2003.
  * [38] W. Ma, A. Trusina, H. El-Samad, W. A. Lim, and C. Tang, “Defining network topologies that can achieve biochemical adaptation.” _Cell_ , vol. 138, no. 4, pp. 760–773, 2009.
  * [39] R. P. Araujo and L. A. Liotta, “Design principles underlying robust perfect adaptation of complex biochemical networks,” _Methods in Molecular Biology_ , (in press).
  * [40] C. P. Robert, _The Bayesian Choice_. Springer Verlag, 2007.
  * [41] P. D. W. Kirk, A. C. Babtie, and M. P. H. Stumpf, “Systems biology (un)certainties.” _Science_ , vol. 350, no. 6259, pp. 386–388, 2015.
  * [42] P. D. W. Kirk, T. W. Thorne, M. P. H. Stumpf, and M. P. H. Stumpf, “Model selection in systems and synthetic biology.” _Current opinion in biotechnology_ , vol. 24, no. 4, pp. 767–774, 2013.
  * [43] T. Toni, D. Welch, N. Strelkowa, A. Ipsen, M. P. H. Stumpf, and M. P. H. Stumpf, “Approximate bayesian computation scheme for parameter inference and model selection in dynamical systems.” _Journal of the Royal Society Interface_ , vol. 6, no. 31, pp. 187–202, 2009.
  * [44] R. D. Wilkinson, “Approximate bayesian computation (abc) gives exact results under the assumption of model error.” _Statistical applications in genetics and molecular biology_ , vol. 12, no. 2, pp. 129–141, 2013.
  * [45] C. P. Robert, J.-M. Cornuet, J.-M. Cornuet, J.-M. Marin, J.-M. Marin, N. S. Pillai, and N. S. Pillai, “Lack of confidence in approximate bayesian computation model choice.” _Proceedings of the National Academy of Sciences of the United States of America_ , vol. 108, no. 37, pp. 15 112–15 117, 2011.
  * [46] C. P. Barnes, T. W. Thorne, S. Filippi, and M. P. H. Stumpf, “Considerate approaches to constructing summary statistics for abc model selection,” _Statistics and Computing_ , vol. 22, no. 6, pp. 1181–1197, 2012.
  * [47] D. Prangle, P. Fearnhead, M. P. Cox, P. J. Biggs, and N. P. French, “Semi-automatic selection of summary statistics for abc model choice.” _Statistical applications in genetics and molecular biology_ , vol. 13, no. 1, pp. 67–82, 2014.
  * [48] E. J. Crampin, E. A. Gaffney, and P. K. Maini, “Reaction and diffusion on growing domains: scenarios for robust pattern formation,” _Bulletin of Mathematical Biology_ , vol. 61, pp. 1093–1120, 1999.
  * [49] J. W. Little, “Robustness of a gene regulatory circuit,” _The EMBO Journal_ , vol. 18, pp. 4299–4307, 1999.
  * [50] L. Ma and P. A. Iglesias, “Quantifying robustness of biochemical network models,” _BMC Bioinformatics_ , vol. 3, p. 38, 2002.
  * [51] M. Morohashi, A. E. Winn, M. T. Borisuk, H. Bolouri, J. Doyle, and H. Kitano, “Robustness as a measure of plausibility in models of biochemical networks,” _Journal of Theoretical Biology_ , vol. 216, pp. 19–30, 2002\.
  * [52] H. Kitano, “Cancer as a robust system: implications for anticancer therapy,” _Nature Reviews Cancer_ , vol. 4, pp. 227–235, 2004.
  * [53] T. Tian, “Robustness of mathematical models for biological systems,” in _Proceedings of the 11th Computational Techniques and Applications Conference (CTAC-2003)_ , J. Crawford and A. J. Roberts, Eds., vol. 45, 2004, pp. C565–C577.
  * [54] D. A. Rand, “Mapping global sensitivity of cellular network dynamics: sensitivity heat maps and a global summation law,” _Journal of The Royal Society Interface_ , vol. 5, pp. S59–S69, 2008.
  * [55] J. Raerinne, “Robustness and sensitivity of biological models,” _Philosophical Studies_ , vol. 166, pp. 285–303, 2013.
  * [56] L. Li, M. Han, and Y. Liu, “Existence and uniqueness of traveling wave front of a nonlinear singularly perturbed system of reaction-diffusion equations with a Heaviside step function,” _Journal of Mathematical Analysis and Applications_ , vol. 410, pp. 202–212, 2014.
  * [57] Y. Matsuoka, N. Jahan, and H. Kurata, “S-system-based analysis of the robust properties common to many biochemical network models,” _Bioprocess and Biosystems Engineering_ , vol. 39, pp. 735–746, 2016.
  * [58] J. Stelling, U. Sauer, Z. Szallasi, F. J. Doyle, and J. Doyle, “Robustness of cellular functions,” _Cell_ , vol. 118, pp. 675–685, 2004.
  * [59] F. Blanchini and E. Franco, “Structurally robust biological networks,” _BMC Systems Biology_ , vol. 5, p. 74, 2011.
  * [60] F. Tostevin, P. R. ten Wolde, and M. Howard, “Fundamental limits to position determination by concentration gradients,” _PLoS Computational Biology_ , vol. 3, p. e78, 2007.
  * [61] A. Wagner, “Robustness and evolvability: a paradox resolved,” _Proceedings of the Royal Society B: Biological Sciences_ , vol. 275, pp. 91–100, 2008.
  * [62] J. T. Young, T. S. Hatakeyama, and K. Kaneko, “Dynamics robustness of cascading systems,” _PLOS Computational Biology_ , vol. 13, p. e1005434, 2017.
  * [63] M. Ó Searcóid, _Metric Spaces_. Springer London, 2007.
  * [64] S. Chaudhuri, S. Gulwani, and R. Lublinerman, “Continuity and robustness of programs,” _Communications of the ACM_ , vol. 55, pp. 107–115, 2012.
  * [65] R. Majumdar, E. Render, and P. Tabuada, “Robust discrete synthesis against unspecified disturbances,” in _Proceedings of the 14th International Conference on Hybrid Systems: Computation and Control_ , ser. HSCC ’11. Association for Computing Machinery, 2011, pp. 211–220.
  * [66] G. R. Wood and B. P. Zhang, “Estimation of the Lipschitz constant of a function,” _Journal of Global Optimization_ , vol. 8, pp. 91–103, 1996.
  * [67] J. Zabczyk, _Mathematical Control Theory_ , 2nd ed. Birkhäuser, 2020.
  * [68] M. Secrier, T. Toni, M. P. H. Stumpf, and M. P. H. Stumpf, “The abc of reverse engineering biological signalling systems.” _Molecular Biosystems_ , vol. 5, no. 12, pp. 1925–1935, 2009.
  * [69] A. C. Babtie, P. D. W. Kirk, and M. P. H. Stumpf, “Topological sensitivity analysis for systems biology.” _Proceedings of the National Academy of Sciences of the United States of America_ , vol. 111, no. 52, pp. 18 507–18 512, 2014.
  * [70] T. Toni, M. P. H. Stumpf, and M. P. H. Stumpf, “Simulation-based model selection for dynamical systems in systems and population biology.” _Bioinformatics (Oxford, England)_ , vol. 26, no. 1, pp. 104–110, 2010.
  * [71] C. P. Barnes, D. Silk, M. P. H. Stumpf, and M. P. H. Stumpf, “Bayesian design strategies for synthetic biology,” _Interface Focus_ , vol. 1, no. 6, pp. 895–908, 2011.
  * [72] Y. Schaerli, A. Munteanu, M. Gili, J. Cotterell, J. Sharpe, and M. Isalan, “A unified design space of synthetic stripe-forming networks,” _Nature communications_ , vol. 5, p. 4905, 2014.
  * [73] O. Purcell, B. Jain, J. R. Karr, M. W. Covert, and T. K. Lu, “Towards a whole-cell modeling approach for synthetic biology,” _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , vol. 23, no. 2, p. 025112, 2013\.