## Abstract

Biomolecular circuits with two distinct and stable steady states have been identified as essential components in a wide range of biological networks, with a variety of mechanisms and topologies giving rise to their important bistable property. Understanding the differences between circuit implementations is an important question, particularly for the synthetic biologist faced with determining which bistable circuit design out of many is best for their specific application. In this work we explore the applicability of Sturm's theorem—a tool from nineteenth-century real algebraic geometry—to comparing ‘functionally equivalent’ bistable circuits without the need for numerical simulation. We first consider two genetic toggle variants and two different positive feedback circuits, and show how specific topological properties present in each type of circuit can serve to increase the size of the regions of parameter space in which they function as switches. We then demonstrate that a single competitive monomeric activator added to a purely monomeric (and otherwise monostable) mutual repressor circuit is sufficient for bistability. Finally, we compare our approach with the Routh–Hurwitz method and derive consistent, yet more powerful, parametric conditions. The predictive power and ease of use of Sturm's theorem demonstrated in this work suggest that algebraic geometric techniques may be underused in biomolecular circuit analysis.

## 1. Introduction

The field of synthetic biology has rapidly matured to the point where it is now possible to produce complex synthetic networks with prescribed functions and levels of performance [1]. As in other fields of engineering, advances have been enabled by the use of small interchangeable modules that are ‘functionally equivalent’ from an input–output perspective [2]. Bistable circuits—which play a role in essential biological processes, including cell fate specification [3], cell cycle progression [4] and apoptosis [5]—make up a particularly large and diverse functionally equivalent set [6]. Effectively characterizing and comparing these biocircuits is crucial for determining which design is in some sense optimal for a particular context.

Ordinary differential equation (ODE) models can be powerful tools for identifying and contrasting biocircuits' ‘dynamic phenotypes' [7]. Numerical simulations are often used in the analysis of these models; however, analytical criteria that focus on topology can provide a more exact assessment of a circuit's properties [8,9]. (Here, we use *topology* to mean the particular set of interactions between regulatory parts.) A novel analytical tool that can provide topology-based insights can be found in Sturm's theorem [10], developed in 1835 as a solution to the problem of finding the number of real roots of an algebraic equation with real coefficients over a given interval. Despite its predictive power, this ‘gem of 19th century algebra and one of the greatest discoveries in the theory of polynomials' [11] remains an unexploited tool for analysis of biological circuit models.

In this work, we demonstrate an approach to bistable circuit discrimination based on Sturm's theorem that can give the boundaries of the regions of bistability as exact analytic expressions, eliminating the need for numerical simulation. We compare the regions of bistability for two variants of the classic double-negative toggle switch as well as two positive feedback circuits, one of which is based on the bacteriophage *λ* promoter P_{RM}. We then show that while a purely monomeric version of the genetic toggle cannot be bistable, a single competitive activating species added to the circuit leads to bistability in a non-contiguous region of parameter space. Lastly, we use a model of an RNA aptamer-based bistable switch to compare our Sturm's theorem approach to another based on the control theoretic Routh–Hurwitz method. Overall, our results highlight a new use for Sturm's theorem for identifying potential differences between functionally equivalent bistable biocircuits, and serve as a (re-)introduction to the method as a general tool for studying the kinds of polynomial expressions that often arise when modelling biological systems.

## 2. Mathematical preliminaries

Sturm's theorem gives the number of distinct real roots of a univariate polynomial *f*(*x*) in a particular interval. To apply the theorem, we must first construct the *Sturm sequence*, a set of polynomials ℱ = {*f*_{0}, *f*_{1}, …, *f*_{m}} defined as
where rem(*f*_{i}, *f*_{i+1}) is the remainder of the polynomial long division of *f _{i}* by

*f*

_{i+ 1}. The sequence ends at

*f*, when

_{m}*f*

_{m− 1}divided by

*f*gives a remainder of zero. For a polynomial of degree

_{m}*n*, there are

*m*≤

*n*+ 1 Sturm polynomials in the sequence.

### Theorem 2.1

(**Sturm's theorem**). *Let f*(*x*) *be a real-valued univariate polynomial and* *a*, *b* ∈ ℝ ∪ {−∞, +∞}, *with a* < *b and* *f*(*a*), *f*(*b*) ≠ 0. *Then the number of zeroes of f*(*x*) *in the interval* (*a, b*) *is the difference*
*where* ℱ *is the Sturm sequence of f*(*x*), *and the* variations var(ℱ, *a*) *and* var(ℱ, *b*) *are the number of times that consecutive non-zero elements of the Sturm sequence—evaluated at a and b, respectively—have opposite signs.* (*Adapted from* [12]*.*)

Our approach involves identifying regions of bistability by finding conditions that lead to three steady states, without requiring numerical determination of the exact values or stability of the equilibrium points. While it is, in general, not possible to draw conclusions on the stability properties of equilibria by simply counting their number, the circuits under consideration enjoy important properties—namely they are dissipative, their linearizations are positive and the characteristic polynomials of their Jacobians have all non-constant terms positive—that allow us to relate their degree and number of equilibria to the stability properties of each equilibrium (see the electronic supplementary material). For such circuits, when three equilibria are present, two of them must be stable and one must be unstable.

## 3. Overview of the method

To determine the region in parameter space in which a particular circuit exhibits bistability, we first find the polynomial *f*(*x*) that describes the equilibrium state of the system (where *x* represents the concentration of a circuit species of interest) and construct its Sturm sequence ℱ = {*f*_{0}, *f*_{1}, …, *f*_{m}} as described above. Because we are interested in all possible positive concentrations, *x* → 0 and *x* → +∞ are chosen as the limits of the interval on which Sturm's theorem is applied. For each limit, we generate different sets of polynomial inequalities by independently setting each element of ℱ to be either greater than or less than zero; however, only some of these inequality sets may be combined in such a way to yield a variation difference var(ℱ, 0) − var(ℱ, +∞) = 3 and thus three steady states. (For example, with var(ℱ, 0) = 4, only those complementary sets with var(ℱ, +∞) = 1 need to be considered.) These combinations of inequalities are enumerated and then tested for *logical consistency*, i.e. whether all inequalities can be simultaneously satisfied. It is important to note that, because Sturm's theorem is only concerned with consecutive non-zero elements and zeroes are ignored, a particular Sturm polynomial may be equal to zero without affecting the total number of sign changes. In cases where zeroes are valid options in the sequence, strict inequalities are made non-strict (e.g. ‘>’ → ‘≥’).

Symbolic manipulation is done using Mathematica; a sample Mathematica notebook in which the method is applied is included as electronic supplementary material.

## 4. Results

### 4.1. Genetic toggle circuits

A recent study identified a set of 11 *minimal bistable networks* (MBNs), simple two-gene circuits with the capacity for bistability that do not also contain a smaller bistable subcircuit [13]. One of these MBNs, a double-negative toggle switch consisting of two dimeric repressors (figure 1*a*, top), was among the very first synthetic biocircuits built and modelled [14]. This *dimer–dimer* (DD) toggle design has since gone on to be used in a wide range of synthetic biological applications, including the manipulation of fluxes of the *Escherichia coli* metabolic network [15] and as part of a circuit involved in programmed autonomous cellular diversification [16]. A second MBN of particular interest, herein referred to as the *monomer–dimer* (MD) toggle, is a double-negative switch variant in which one of the repressors functions as a monomer (figure 1*a*, bottom). While to the best of our knowledge no MD toggle circuit has been constructed, the components necessary for its implementation exist in the form of monomeric single-chain transcriptional repressor proteins [17], as well as in transcription-activator-like effectors (TALEs) and CRISPR/Cas nucleases that have been engineered as repressors [18,19]. Other exotic toggle-like circuit topologies have also been proposed and/or built [13,18,20–22]. An understanding of the differences in how these various toggles perform can be beneficial, in particular for those circuits that have yet to be built and for which *a priori* knowledge of their behaviour could aid in their development.

As a first demonstration of our approach, we apply Sturm's theorem to the DD and MD toggle circuits. Beginning with a chemical reaction network formulation and assuming mass-action kinetics we derive ODE models of the two toggles (electronic supplementary material, equations (S1) and (S2)). At equilibrium, the concentrations of *P*_{1} and *P*_{2} in the MD system are given by
4.1and in the DD system,
4.2where *X*_{itot} is the total concentration of gene *i*, *β*_{i} = *k*_{basi}/*k*_{degi} is the ratio of basal production and degradation rate constants for protein *i*, *K*_{xd} = {*K*_{md}, *K*_{dd}} is the Michaelis constant for *P*_{1}, and *K*_{2} is the Michaelis constant for *P*_{2}. (Recall that Michaelis constants represent the repressor concentrations that yield 50% of the maximum production rate of their target proteins and are determined by the protein–DNA and dimerization dissociation constants; see electronic supplementary material for details.) Systems (4.1) and (4.2) may be written in terms of *P*_{1eq} alone, as
4.3and
4.4With the following scaling of the DNA and protein concentrations:
4.5and
4.6we may write equations (4.3) and (4.4) as non-dimensional polynomials in :
4.7and
4.8for the MD and DD toggles, respectively. Every positive root of these equilibrium polynomials also gives a positive steady-state concentration for every other circuit component.

We may now determine the regions of bistability in the plane of total (non-dimensionalized) DNA concentrations and . As outlined in §3, the first step is the construction of a Sturm sequence ℱ for each equilibrium polynomial. These sequences are given in the electronic supplementary material. It can be seen that, for the DD toggle, some of the polynomials in the sequence contain a term that renders them indeterminate when . (The same term is responsible for one of the polynomials going to zero, thus terminating the sequence prematurely.) To avoid these problems, it is necessary to generate two different sequences for the DD toggle: one for which we assume in equation (4.8), and another with . In this latter case, the DD toggle equilibrium polynomial equation (4.8) becomes 4.9

Following construction of the toggle Sturm sequences, we evaluate each at and and find the conditions that give a variation difference of 3. The MD toggle Sturm sequence ℱ_{md} has a maximum possible variation of 3 and only one combination of inequalities that can give rise to bistability: when var(ℱ_{md}, 0) = 3 and var(ℱ_{md}, +∞) = 0. In contrast, the DD toggle sequence ℱ_{dd} could in principle yield as many as five positive steady states; however, only three are admitted as there are no combinations of inequalities that have a variation difference of 5 or 4 and are logically consistent. (See §3 for details.) Inequality sets for the DD toggle are listed (in the compact form {± ± · · · ±}) in tables 1 and 2.

The valid inequality sets were reduced with Mathematica to yield analytic expressions for the regions of bistability in the plane of and . For the MD toggle, this bistability region is found to be 4.10with . The DD toggle inequality sets combine to give a continuous region of bistability that is most easily written as the intersection of 4.11and 4.12It is important to note that, with the exception of the DD toggle at (which was treated by analysing a second Sturm sequence, as described above), none of the potential zeroes in the Sturm polynomial denominators required special treatment nor did they present any problems in determining the regions of bistability.

Our analysis shows that the DD toggle operates as a bistable switch over a substantially greater range of (non-dimensionalized) DNA concentrations than does the MD toggle (figure 1*b*), indicating that the DD topology is more functionally robust to variations in DNA concentrations and rate parameters. Furthermore, the DD switch can operate with significantly lower values of and : > 75% and > 50% less, respectively. Worth highlighting is the fact that, through application of our method to the non-dimensionalized versions of the original systems, we have demonstrated that the *function* of bistability is determined by a particular combination of the system parameters—gene concentrations, protein basal production and degradation rate constants, and Michaelis constants—and thus there is no single parameter more important than any other for achieving bistability.

### 4.2. Computational support

The approach presented here does not require any numerical simulation; however, some computational validation of our results may be considered of value. For both toggle circuits, and for each of the valid combinations of sign(ℱ, 0) and sign(ℱ, +∞), 1000 random values of and were selected from inside and outside of the predicted bistable regions and plugged in to the appropriate equilibrium polynomial (equations (4.7) or (4.8)) that were then solved numerically. In all cases, the number of equilibria found matched the number determined by Sturm's theorem: three equilibria were found inside the bistable regions and only one equilibrium was found outside.

We may also check the stability of the various steady states using the circuits' Jacobian *J* and characteristic polynomial *p _{J}*(

*λ*) = det(

*λ*−

**I***J*). It was recently shown that if all off-diagonal components of the Jacobian are non-negative (i.e. it is a Metzler matrix), or if the Jacobian may be transformed to have such a form, then any equilibrium is unstable if and only if the constant term of

*p*(

_{J}*λ*) has a sign opposite to that of all other terms in

*p*(

_{J}*λ*) [20]. We use this condition on the constant term of

*p*(

_{J}*λ*) to confirm that each bistable solution set contains one and only one unstable steady state.

The inequalities that satisfy the *p _{J}*(

*λ*) constant term condition are 4.13and 4.14for the MD and DD toggles, respectively. For each bistable solution set found, we substituted the values and into equations (4.13) and (4.14) and confirmed that only one of the three solutions satisfies the appropriate instability condition.

It is worth noting that the time required to test pairs scales linearly with the number of pairs, so while testing a small number can be done relatively quickly, as the number of pairs becomes appreciable the time can be significant—up to 3 h to test 600 000 random values of and .

### 4.3. Bistable single-gene circuits

The single gene system consisting of bacteriophage *λ* repressor and its promoter P_{RM} (with its three operator sites OR1, OR2 and OR3) also exhibits bistability [23]. Using the dimensionless model given in [23], we have that the steady state concentration of protein satisfies
4.15where *γ* is the rescaled degradation rate constant, *α* represents the increase in protein production resulting from dimer binding to OR2, and *σ*_{1} and *σ*_{2} are the relative (to OR1) affinities for OR2 and the negatively regulating OR3, respectively. (For simplicity, we set the gene copy number equal to one.) With *σ*_{1} = 2 and *σ*_{2} = 0.08 [23], the associated Sturm sequence has only two inequality sets with and one set with that are logically consistent and together give bistability (table 3).

We can compare the bistability region of the multi-operator P_{RM} system with that of a simple positive feedback circuit consisting of a dimeric repressor and only one operator site (MBN *kqw* in [13]). Rescaled as in equation (4.15), we have
4.16

As with the MD toggle, this polynomial also has a maximum possible variation of 3 and thus only a single combination of inequalities that give rise to bistability, in the region given by 4.17

The bistable regions (in *α* – *γ* space) for these single gene circuits are shown in figure 2*a*. (The one-dimensional regions of bistability shown in bifurcation diagrams in [23,24] are plotted for comparison.) It can be seen that the *λ* repressor circuit is bistable over a larger range and with lower values of the degradation rate constant. Interestingly, with *α* = 11 ([23] and references therein), the single operator circuit would just barely function as a bistable circuit, and any small fluctuation in circuit parameters would render it non-functional.

We can also use our method to determine how the strength of the negative feedback (*σ*_{2}) affects bistability. Keeping *σ*_{1} = 2, and using *α* = 11 and *γ* = 4.5 (centred in the bistable region at *α* = 11; figure 2*a*), we find that *σ*_{2} can increase 12-fold to approximately 0.96 before bistability is lost. In general, significant increases in *σ*_{2} require similar increases in *α* for bistability to be maintained, with the range of allowable *γ* narrowing as a result (figure 2*b*).

### 4.4. Purely monomeric toggle circuits

In a recent study [18], it was shown that a toggle variant consisting of two monomeric repressors derived from TALE protein DNA-binding domains cannot support bistability, but that the introduction of positive feedback in the form of two TALE-based activators that compete with the repressors for promoter access makes the system bistable. However, as we show below, it is not necessary to have two symmetric positive feedback loops for bistability; a single activator added to the monomeric toggle is sufficient.

We begin by combining components from the ‘mutual repressor’ and ‘competitive feedback’ models of Lebar *et al*. [18], so that the two types of promoters—one containing a single binding site for both an activator and a repressor, and a second that is only responsive to a repressor—are present in one hybrid circuit:
4.18

In keeping with the nomenclature of [18], [*T*_{AV}] represents the concentration of the activator species, and [*T*_{AK}] and [*T*_{BK}] are the concentrations of the two repressor species. (For simplicity, we have not included terms related to the fluorescent reporters as they have no effect on the capacity for bistability.) So that our results may be easily comparable to those in [18], we include a basal rate constant *k*_{b} to account for expression from the pristinamycin- and erythromycin-inducible promoters incorporated in the original system for external control. Other parameters of the models in [18] set to either zero or one in their simulations are, for clarity of presentation, also set to those same values here.

Equation (4.18) may be combined into a single equilibrium polynomial in [*T*_{AK}]:
4.19

The maximum possible variation that the Sturm sequence associated with equation (4.19) may have is 3 and there is thus only a single combination of inequalities that can yield three steady states.

To establish the relationship between the repressor and activator expression strengths (*k _{K}* and

*k*, respectively) and the capacity for bistability, we set

_{V}*k*

_{b}to the uninduced steady-state value of approximately 0.04 M h

^{−1}and the degradation rate constant

*k*

_{d}= 0.1 h

^{−1}(values taken from [18]). Interestingly, we find two unconnected bistable subregions: one in which the repressor promoter strength must be greater than 10 times that of the activator, and another in which the activator promoter strength must be greater than approximately 2 times that of the repressor (figure 3). Thus, while a purely monomeric toggle without positive feedback cannot be bistable—indeed, the equilibrium polynomial for such a circuit is only second order and could never have more than two real roots—the addition of a single activator with competitive feedback is all that is needed for bistability. Also noteworthy is that the values of

*k*and

_{K}*k*that give bistability in the simulation of the full ‘competitive feedback’ model of [18] (180 and 30 h

_{V}^{−1}, respectively) only give a monostable response in this hybrid circuit.

### 4.5. Sturm's theorem and the Routh–Hurwitz stability criterion: a comparison of analytical methods

The Routh–Hurwitz stability criterion is commonly employed by control theorists to test if a polynomial admits positive roots. Typically applied to characteristic polynomials of linear time invariant systems to establish their stability, the method can also be used to find the parameter range in which a polynomial equilibrium condition admits a desired number of positive roots. Indeed, Sturm's theorem and the Routh–Hurwitz method are related, and the validity of the latter can be demonstrated using Sturm chains. It may be argued that building a Routh table is a simpler procedure than applying Sturm's theorem because it only requires basic arithmetic operations on the polynomial coefficients; however, its relative simplicity comes at a price: information about the real or complex nature of the positive roots is lost. In contrast, Sturm's theorem allows for an exact determination of the number of real positive roots, as well as whether or not they are physically meaningful for a biochemical system.

A recent study of a novel biomolecular circuit that uses modulation of enzyme activity by inhibitory RNA aptamers to achieve bistability [20] provides us with an opportunity to compare our approach with one based on the Routh–Hurwitz method. Interestingly, like the monomeric circuit from [18] and its asymmetric variant described above, and proposed MBN *bcdh* in [13], this aptamer-based circuit achieves bistability with monomeric regulation.

The equilibrium condition for the concentration of one of the enzyme species in its active form (given here by *x*; see electronic supplementary material) is a complicated polynomial of fourth order:
where *E*^{tot} is the total concentration of each of the two enzymes used in the circuit and Greek letters indicate reaction rates.

The inequality sets for this circuit are listed along with their allowabilities in table 4. Combining these inequalities, we find that for bistability it is necessary that 4.21

These conditions are identical to those identified in [20] by applying the Routh–Hurwitz table together with imposing conditions to guarantee *real* positive roots and conditions on the constant term of the characteristic polynomial to achieve bistability. Additional details are provided in the electronic supplementary material.

## 5. Discussion

In this work, we have for the first time applied Sturm's theorem to the study of biological circuits; specifically, biocircuits that can adopt two stable steady states. The primary results of our analyses are exact analytical expressions defining the various circuits' regions of bistability, and without the need for parameter surveys and numerical simulation. We take advantage of the fact that the bistable circuits considered here are dissipative, positive in their linearization, and have characteristic polynomials with all non-constant terms positive, properties that allow us to ascertain the stability of their equilibria without computation. It is important to emphasize that these are not unique features of these particular circuits and that many real systems share these properties. For any such system, once the number of steady-states is determined—by Sturm's theorem or by other methods—the stabilities are also known.

So following our analysis, can it be said that some topologies are in some sense ‘better’ than others? Between the two different genetic toggle variants, we see that when rate parameters are fixed, the toggle consisting of two dimeric repressor species functions as a bistable switch over a wider range of DNA concentrations than one composed of one dimeric and one monomeric repressor. This result provides a strong motivation for choosing a DD toggle over a MD toggle in any application where there is considerable uncertainty or variability in parameter values or DNA concentrations (e.g. when DNA is in the form of plasmids without strict copy-number control). In single-gene positive feedback circuits exhibiting bistability, we see benefits to having additional operator sites: without both OR1 and OR2 in the *λ* repressor–P_{RM} system, the enhancement *α* = 11 would barely be sufficient for bistability. Interestingly, although feedback at OR3 is negative, it is not strong enough to significantly affect the circuit's ability to function as a bistable switch. Taken together, these two results suggest that the promoter architecture of the *λ* system may have evolved to allow for both robust bistability due to the positive feedback as well as reduced variability or other benefit of the small negative autoregulation.

Although relatively little known within the biological sciences, Sturm's theorem has found applicability in a number of other areas where polynomials play an important role, including computational mathematics [25], dynamical systems [26,27], robotics [28] and finance [29]. Additional biological applications are possible, for example, as a tool to predict a large number of new bistable topologies or rule out those that do not have the capacity for bistability (such as chemical reaction network theory, previously [13,24]).

Other methods from algebraic geometry, including the closely related Descartes' rule of signs, have also found use in biological systems analysis, with recent applications to model discrimination [30] and the study of chemical reaction networks [31,32] and multisite phosphorylation systems [33], among many others. (This latter paper is particularly relevant to our own work, as it also highlights the benefits of treating parameters symbolically rather than numerically.) For the specific task required in this work—determining the number of positive real roots of a biocircuit's equilibrium polynomial—Descartes' rule of signs may be simpler to apply than Sturm's theorem; however, it is definitively less powerful in that it can only give an upper bound on the number of real roots.

We note that not all bistable circuits of interest may be easily studied using our approach. First, it may not be possible to capture the circuit's equilibrium state with a single polynomial. When an equilibrium polynomial can be found, those of third order are the simplest to study, because there is only one combination of inequalities that can give rise to bistability. As the order increases, the number of inequality sets to test and the complexity of the Sturm sequences can also increase to the point where application of the theorem is impractical. However, this is highly dependent on the complexity of the polynomial coefficients and the number of Sturm sequences of fixed sign. Indeed, as we demonstrated with our analysis of the *λ* repressor–P_{RM} system, even a seventh-order polynomial can be tractable. For particularly large circuits for which an equilibrium polynomial cannot be derived, or for which the Sturm sequence is overly complex, a decomposition of the system into smaller modules can be beneficial.

The analysis of an equilibrium polynomial with a zero root—for example, that describing a recently developed DNA-based circuit that achieves bistability through a combination of mutual inhibition and autocatalysis [34]—presents challenges, because Sturm's theorem requires that neither of the limits of the region of interest (in our case, 0 and +∞) be roots. Under these circumstances, one could choose an *ε* < 0 for the lower limit; however, in this case, there is no guarantee that the three roots found by Sturm's theorem will all represent physically real, positive quantities (because one or more of the roots could in theory lie between *ε* and 0).

### 5.1. Extensions to Sturm's theorem and related methods

The Sturm's theorem approach described here can be directly applied to systems whose steady state may be characterized by a univariate polynomial with integer exponents and simple roots. While such systems are common, examples of equilibrium polynomials not of this particular form can also be found. For example, exponents that are typically given integer values in simplified models and/or to reflect strong multimerization (as in this work) or other reactions with significant positive cooperativity [35] can be found to have fractional values when the models are fitted to empirical data [36–38]. In limited cases, these generalized polynomials can be turned to a proper polynomial with a simple substitution that does not affect the number of zeros (e.g. *u* = *x*^{1/N}, if all exponents of *x* are multiples of 1/*N*). There are also simple extensions to Sturm's theorem that deal with multiple roots [39] and multivariate models [12,40]. Should only an upper bound on the number of real roots of multivariate and generalized polynomial models be desired, extensions to the related Descartes' rule of signs are also available [41–43].

## Authors' contributions

The project was conceived by D.S.-G., E.F. and R.M.M. D.S.-G. and E.F. performed the mathematical analyses, with computational support provided by T.Z. The paper was written by D.S.-G. and E.F. and edited by all of the coauthors.

## Competing interests

We declare we have no competing interests.

## Funding

This research is funded in part by the National Science Foundation through grant CMMI 1266402, and the Gordon and Betty Moore Foundation through grant no. GBMF2809 to the Caltech Programmable Molecular Technology Initiative.

## Acknowledgements

A large number of people contributed to this work with insights and comments. The authors particularly thank Andras Gyorgy, Yutaka Hori, Scott C. Livingston, Anne Shiu, Eduardo Sontag, Elisenda Feliu, Zvi H. Rosen, Jaap Top and Brian Ingalls.

- Received March 31, 2015.
- Accepted June 1, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.