Cooperativity, Absolute Interaction, and Algebraic Optimization

TL;DR

This paper introduces a measure of cooperativity based on minimal absolute interaction, formulates it as an algebraic optimization problem, and applies it to analyze hemoglobin binding behaviors under various conditions.

Contribution

It presents a novel algebraic optimization approach to quantify cooperativity through minimal absolute interaction, using nonlinear algebra tools like SCIP.

Findings

Minimal absolute interactions align with chemical modifications' expected outcomes.

The measure ranks cooperativity differently than traditional Hill slope analysis.

Application to hemoglobins demonstrates the method's practical utility.

Abstract

We consider a measure of cooperativity based on the minimal absolute interaction required to generate an observed titration behavior. We describe the corresponding algebraic optimization problem and show how it can be solved using the nonlinear algebra tool \texttt{SCIP}. Moreover, we compute the minimal absolute interactions for various binding polynomials that describe the oxygen binding of various hemoglobins under different conditions. While calculated minimal absolute interactions are consistent with the expected outcome of the chemical modifications, it ranks the cooperativity of the molecules differently than the maximal Hill slope.

Equations

$$\mathrm{W}=(w_I{)}_{I\subseteq [n]}\in ({\mathbb{R}}_{>0}{)}^{2^n}.$$

$$\Phi (\mathrm{W}):=a_n{\Lambda }^n+\cdots +a_1\Lambda +a_0,$$

$$a_k:=\sum _{|I|=k}\prod _{I^{\prime }\subseteq I}{w}_{I^{\prime }}\in {\mathbb{R}}_{>0}.$$

$$a_0$$

$$a_1$$

$$a_2$$

$$a_3$$

$$\|{\rm W}\|:=\prod_{|I|>1}\max\Big{(}w_{I},w_{I}^{-1}\Big{)}.$$

$$\|P\|:=\min\Big{\{}\|{\rm W}\|\left.\vphantom{\Big{\{}}\;|\;\vphantom{\Big{\}}}\right.{\rm W}\in(\mathbb{R}_{>0})^{2^{n}}\text{ with }\Phi({\rm W})=P\Big{\}}.$$

$$P_1:=$$

$$P_2:=$$

$$P_3:=$$

$$\Vert P_1\Vert =13.50,\Vert P_2\Vert =3,\Vert P_3\Vert =1.$$

$$\begin{array}[]{lll}\displaystyle\operatorname*{minimize}_{{\rm W}}&&\|{\rm W}\|\\[8.53581pt] \operatorname*{subject~{}to}&&\Phi({\rm W})=P\end{array}$$

$$\begin{array}[]{lcl}\displaystyle\operatorname*{minimize}_{{\rm W}}&&\displaystyle\prod_{|I|>1}\max\Big{(}w_{I},w_{I}^{-1}\Big{)}\\[8.53581pt] \operatorname*{subject~{}to}&&a_{1}=w_{1}+w_{2}+w_{3}+w_{4},\\[2.84526pt] &&a_{2}=w_{1}w_{2}w_{12}+w_{1}w_{3}w_{13}+w_{1}w_{4}w_{14}\\ &&\qquad\qquad+w_{2}w_{3}w_{23}+w_{2}w_{4}w_{24}+w_{3}w_{4}w_{34},\\[2.84526pt] &&a_{3}=w_{1}w_{2}w_{3}w_{12}w_{13}w_{23}w_{123}+w_{1}w_{2}w_{4}w_{12}w_{14}w_{24}w_{124}\\ &&\qquad\qquad+w_{1}w_{3}w_{4}w_{13}w_{14}w_{34}w_{134}+w_{2}w_{3}w_{4}w_{23}w_{24}w_{34}w_{234},\\[2.84526pt] &&a_{4}=w_{1}w_{2}w_{3}w_{4}w_{12}w_{13}w_{14}w_{23}w_{24}w_{34}w_{123}w_{124}w_{134}w_{234}w_{1234}.\end{array}$$

$$O_{\mathcal{I}}:=\big{\{}(w_{I})_{I\subseteq[n]}\mid\forall|I|>1:w_{I}\geq 1\;\forall I\notin\mathcal{I}\text{ and }w_{I}\leq 1\;\forall I\in\mathcal{I}\big{\}}$$

$$({\mathbb{R}}_{>0}{)}^{2^n}=\bigcup _{\mathcal{I}\subseteq \{I\subseteq [n]\mid |I|>1\}}{O}_{\mathcal{I}}.$$

$$\Vert W\Vert =\prod _{I\in \mathcal{I}}{w}_I^{-1}\prod _{\begin{array}{c}|I|>1\\ I\notin \mathcal{I}\end{array}}{w}_I=:f_{\mathcal{I}}(W)\text{for }\mathrm{W}\in {\mathcal{O}}_{\mathcal{I}}.$$

$$\begin{array}[]{lll}\displaystyle\operatorname*{minimize}_{{\rm W}\in\mathcal{O}_{\mathcal{I}}}&&f_{\mathcal{I}}({\rm W})\\[8.53581pt] \operatorname*{subject~{}to}&&\Phi({\rm W})=P.\end{array}$$

$$S_n\times ({\mathbb{R}}_{>0}{)}^{2^n}⟶({\mathbb{R}}_{>0}{)}^{2^n},(\sigma ,(w_I))⟼\sigma \cdot (w_I):=(w_{\sigma (I)}).$$

$$\min \{\Vert \mathrm{W}\Vert \mid \mathrm{W}\in O_{\{12\}}\text{ with }\Phi (\mathrm{W})=P\}$$

$$=\min \{\Vert \mathrm{W}\Vert \mid \mathrm{W}\in O_{\{ij\}}\text{ with }\Phi (\mathrm{W})=P\}.$$

$$\underset{W}{\mathrm{minimize}}\underset{\mathcal{I}\subseteq \{I\subseteq [n]\mid |I|>1\}}{\max }{f}_{\mathcal{I}}(W)$$

$$\mathrm{subject}\mathrm{to}\Phi (W)=P.$$

$$\underset{W}{\mathrm{minimize}}r$$

$$\mathrm{subject}\mathrm{to}\{\begin{array}{l}{f}_{\mathcal{I}}(W)\le r\text{for all }\mathcal{I}\subseteq \{I\subseteq [n]\mid |I|>1\},\\ \Phi (W)=P.\end{array}$$

$$f_{\mathcal{I}}(W)=w_{123}^{-1}\cdot \prod _{\begin{array}{c}|I|>1,\\ I\ne \{1,2,3\}\end{array}}{w}_I\text{becomes}{f}_{\mathcal{I}}\circ {\phi }^{-1}(S)=\frac{s_{1234}{s}_{12}^2{s}_{13}^2{s}_{23}^2}{s_{123}^2{s}_1^3{s}_2^3{s}_3^3{s}_4}.$$

$$\underset{S}{\mathrm{minimize}}r$$

$$\mathrm{subject}\mathrm{to}\{\begin{array}{l}{f}_{\mathcal{I}}\circ {\phi }^{-1}(S)\le r\text{for all }\mathcal{I}\subseteq \{I\subseteq [n]\mid |I|>1\},\\ {\sum }_{|I|=k}{s}_I=a_k\text{ for all }k=1,\dots ,n.\end{array}$$

$$(1,0,\dots ):s_1\cdot \prod _{i\ne 1}{t}_i,(0,1,0,\dots ):s_2\cdot \prod _{i\ne 2}{t}_i,(1,1,0,\dots ):s_1{s}_2\cdot \prod _{i\ne 1,2}{t}_i,\dots$$

$$S=\sum _{j=1}^k{v}_1^{(j)}\otimes \cdots \otimes v_n^{(j)}$$

$$S=\sum_{j=1}^{k}\Big{(}\prod_{i=1}^{k}t_{i}^{(j)}\Big{)}\cdot\left(\begin{smallmatrix}s_{1}^{(j)}/t_{1}^{(j)}\\ 1\end{smallmatrix}\right)\otimes\dots\otimes\left(\begin{smallmatrix}s_{n}^{(j)}/t_{n}^{(j)}\\ 1\end{smallmatrix}\right),$$

Full text

1. Introduction

Interaction between components is a key concept of biological systems. While characteristics of a system of independent subunits are determined by the description of each element alone, a complex system with interactions is more than the sum of its parts. A classical example of a small biological system with non-trivial interaction is hemoglobin with its four binding sites for oxygen or other ligands such as carbon monoxide [BHK04, Bar13, Hil13]. The overall binding curve, as a function of oxygen concentration at constant temperature (“isotherm”), is much steeper than an overall binding curve generated by a system of independent sites can be. Hemoglobin “prefers” the extreme states of being fully (un)occupied to intermediate states with one, two or three ligands bound. This dependent binding behavior or an “abnormal” steepness of a response curve is described by the concept of “cooperativity”. Cooperativity is not only considered to be important for (oxygen) transport, but for instance also for the formation of multi-protein complexes [RRS17], general signal transduction processes [Sal+16, LFSR09, MSS16] and the regulation of noise [GMD09, MMD13].

Observations on the macroscopic level, that is on how many sites are occupied have often been connected to microscopic behavior that is which sites are occupied and how certain sites interact (for a review see [SL13]). However, it has been shown that in the framework of the grand canonical ensemble of statistical mechanics, different definitions of cooperative binding are not coinciding when general asymmetric systems with more than two binding sites are considered [MDH16, Mar17a]. In particular, macroscopic phenomena which have been attributed to positive cooperativity in literature can be caused by a system being negative cooperative on the microcopic level [MDH16].

For the characterization of the macroscopic binding behavior, roots of the grand partition function (also called the “binding polynomial”) have been the objects of investigations [Bri83, Bri84, Bri85]. In particular, it has been pointed out that the criterion of having non-real roots is generally the relevant aspect indicating that a binding process relies on non-negligible interaction [MDH16]. In case that the binding polynomial has only real roots, these roots define hypothetical single-site-systems whose sum give the initial overall curve [OCU01, OU04, MU13]. The macroscopic behavior of the system could thus theoretically be engineered by using these independent sites which suggests that the potentially present intrinsic interaction is not essential for the macroscopic functionality of the system. Conversely, in case that the binding polynomial has non-real roots, the overall system cannot be represented as a decoupled system consisting of independent “real” binding sites suggesting that the intrinsic interaction is essential for the macroscopic behavior of the system. This qualitative definition raises the question of how to quantify the cooperativity of a system.

The quantification of cooperativity has predominantly been discussed for symmetric systems consisting of physically indistinguishable sites [Abe16, RZYL18, Ron+]. In particular, in the framework of the grand canonical ensemble, the molecule consisting of identical sites can be easily inferred from the coefficients of the binding polynomial and the corresponding interaction energies can be determined [MDH16]. Moreover, it has been proposed to use the absolute minimal interaction which is required to generate the overall binding behavior as a general measure of cooperativity when the binding sites are not considereed to be physically identical [Mar17]. In this general situation without additional side conditions, an observed overall titration behavior can be explained by an infinite number of different hypothetical molecules. A central question to characterize cooperativity can then be what the minimal absolute interaction across all these molecules is and in particular whether it is zero. This quantitative measure of cooperativity has been defined and upper bounds based on polynomial factorization or on the interaction of a symmetric system, have been derived [Mar17]. However, it has not been clear how to explicitly calculate this measure of cooperativity for larger systems.

In this work, we show that looking for the minimal interaction energy required to generate a given overall titration behavior, leads to an algebraic optimization problem which can be tackled using the software SCIP [Gle+18]. We describe the mathematical problem and calculate examples for binding polynomials from native hemoglobins and some modified variants [CRBG86, Ike+83, Ima73]. For our examples, the minimal interaction reduces when the molecule structure of hemoglobin is disturbed. However, our results rank the degree of cooperativity of chemically treated hemoglobines differently than the maximal Hill slope.

2. Recapitulation of the mathematical framework

In this section, we briefly review the model for ligand binding based on the grand canonical ensemble of statistical mechanics and recall the notion of minimal absolute interaction from [Mar17]. We refer the reader to [Ben13, WG90, Hil13a] for detailed exposition on the model.

Definition 2.1** **(Molecule)

A molecule ${\rm W}$ with $n$ sites is a positive real point, whose coordinates are indexed by the subsets of $[n]:=\{1,\ldots,n\}$ and with $w_{\emptyset}=1$:

[TABLE]

We refer to $w_{I}$ as binding energies if $|I|=1$ and interaction energies if $|I|>1$. For the sake of brevity, we abbreviate $w_{\{i_{1},\ldots,i_{r}\}}$ to $w_{i_{1}\ldots i_{r}}$ for $i_{1}<\dots<i_{r}$ (see Fig. 1).

The binding energy $w_{i}$ encodes the difference between free energy of the molecule in the completely unoccupied state and the state with only a ligand bound to site $i$. The interaction energy $w_{i_{1}i_{2}}$ is the difference between free energy when sites $i_{1},i_{2}$ are occupied and the sum of the differences of energy encoded by $w_{i_{1}}$ and $w_{i_{2}}$. Analogously, $w_{i_{1}\dots i_{r}}$ encodes the discrepancy between the energy required for $r$ ligands to bind to the sites $i_{1},\dots,i_{r}$ and the energy already defined by the interactions of lower degree. Note that we may call $w_{i_{1}\dots i_{r}}$ energies even though they actually represent exponentials of the physical quantity “energy” [Ben13, WG90, Sch75].

Definition 2.2** **(Binding Polynomial)

Given a molecule ${\rm W}=(w_{I})_{I\subseteq[n]}$ with $n$ sites, its binding polynomial $\Phi({\rm W})$ is a univariate polynomial of degree $n$,

[TABLE]

whose positive real coefficients $a_{k}$ are given by

[TABLE]

Example 2.3

Let ${\rm W}=(w_{I})_{I\subseteq[3]}\in(\mathbb{R}_{>0})^{2^{3}}$ be a molecule with $3$ sites. The binding polynomial $\Phi({\rm W})=a_{3}\Lambda^{3}+a_{2}\Lambda^{2}+a_{1}\Lambda+a_{0}$ is a real univariate polynomial of degree $3$ whose coefficients are given by

[TABLE]

Note that two different molecules may have the same binding polynomial and thus the same binding curve. Hence, the map $\Phi$ which maps a molecule with $n$ sites to its binding polynomial is not injective.

Definition 2.4** **(Absolute interaction)

The absolute interaction of a molecule ${\rm W}=(w_{I})\in(\mathbb{R}_{>0})^{2^{n}}$ is given by

[TABLE]

Since $w_{I}$ represents the exponential of a difference of actual energies, which can be positive or negative, $\max(w_{I},w_{I}^{-1})$ represents the exponential of the absolute value of that difference. Hence, the absolute interaction $\|{\rm W}\|$ represents the exponential of the sum of the absolute values of all energy differences. The minimal value $\|{\rm W}\|$ can adopt is $1$ which corresponds to a minimal physical interaction energy of [math].

Definition 2.5** **(Minimal absolute interaction)

Given a binding polynomial $P$ of degree $n$, the minimal absolute interaction is

[TABLE]

We call a molecule $W$ minimal, if $\|W\|=\|\Phi(W)\|$.

In other words, the minimal absolute interaction $\|P\|$ is the minimum of all absolute interactions of molecules with a fixed binding polynomial $P$, or equivalently with a fixed binding curve. This notion is well defined by the following proposition.

Proposition 2.6** **([Mar17, §4])

For any univariate polynomial $P$ of degree $n$ with positive coefficients there exists a molecule ${\rm W}\in(\mathbb{R}_{>0})^{2^{n}}$ with $n$ sites such that $\Phi({\rm W})=P$ and $\|{\rm W}\|=\|P\|$.

Since cooperativity is an emerging property of the interaction between the sites and deeply connected to the binding curve, the minimal absolute interaction of a binding polynomial is a natural candidate for quantifying cooperativity. In addition, [Mar17, §4] shows that it has several properties that could be considered desirable in such a quantifier.

Example 2.7

Consider the following binding polynomials:

[TABLE]

A brief computation reveals that the minimal absolute interactions are

[TABLE]

The computation for $\|P_{1}\|$ is explicitly shown in Example 3.2. Figure 2 illustrates the minimal molecules for each binding polynomial. Since $P_{3}$ factorizes into three real linear factors, its minimal molecule has only trivial interactions.

3. The Algebraic Optimization Problem

In this section, we consider the computation of the minimal absolute interaction as an optimization problem with the absolute interaction as objective function and the set of all molecules sharing the same binding polynomial as feasible set:

[TABLE]

This problem seems simple and its concept can be understood easily. However, it becomes quickly complicated if larger systems are considered. For small $n$ it can be worked out by methods of algebraic optimization (e.g. using sums of squares approaches [Las01, Par03] or critical point methods [GS14]). For $n=4$ explicitly, given the binding polynomial $P=a_{4}\Lambda^{4}+a_{3}\Lambda^{3}+a_{2}\Lambda^{2}+a_{1}\Lambda$ we have:

[TABLE]

One intuitive way to compute the minimal absolute interaction $\|P\|$ is to decompose the space of molecules $(\mathbb{R}_{>0})^{2^{n}}$. We look at $2^{2^{n}-n-1}$ regions on which for $|I|>1$ either $w_{I}\geq 1$ or $w_{I}\leq 1$ holds. Explicitly, for each set of interactions $\mathcal{I}\subseteq\{I\subseteq[n]\mid|I|>1\}$, we define a region

[TABLE]

so that

[TABLE]

On each region $O_{\mathcal{I}}$, the absolute interaction $\|W\|$ is a rational function in the interaction energies which we denote by $f_{\mathcal{I}}(W)$,

[TABLE]

Finding the minimal absolute interaction inside $O_{\mathcal{I}}$ is a straight-forward problem in polynomial optimization:

[TABLE]

However, this approach requires solving $2^{2^{n}-n-1}$ polynomial optimization problems, one for each region $\mathcal{O}_{\mathcal{I}}$, which is not easily feasible.

We can reduce the number of regions that require consideration by exploiting the intrinsic symmetry of the problem. To do this, we need to do a small journey in algebra and use the notion of group. In mathematics, a (commutative) group is a set $S$ equipped with a binary operation $\diamond$ that combines any two elements to form a third element of that set and such that this operation is associative (and commutative), the set contains an identity element ($\exists e\in S$ s.t. for all $x\in S$, $e\diamond x=x$) and any element of $S$ is invertible (for all $x\in S$, there exists $y\in S$ s.t. $x\diamond y=e$). The basic example for groups is the set of integers $\mathbb{Z}$ equipped with the binary operation of addition.

In our setting, we want to exploit symmetry. We will consider the symmetric group $S_{n}$ which is the set of all bijections from $\{1,\ldots,n\}$ to $\{1,\ldots,n\}$ equipped with composition of maps. In the sequel, we use the following standard notation: $(ab)$, for $a,b$ in $\{1,\ldots,n\}$, denotes the bijection that permutes $a$ and $b$ in $\{1,\ldots,n\}$.

There exists a natural action of the symmetric group $S_{n}$ on the space of molecules $(\mathbb{R}_{>0})^{2^{n}}$ which permutes the binding sites and the corresponding interaction energies,

[TABLE]

This action permutes the regions $\mathcal{O}_{\mathcal{I}}$ and preserves binding polynomials as well as absolute interactions:

Proposition 3.1

The binding polynomial and the minimal absolute interaction of a molecule are invariant under the group action of $S_{n}$ i.e., $\Phi(\sigma\cdot{\rm W})=\Phi({\rm W})$ and $\|\sigma\cdot{\rm W}\|=\|{\rm W}\|$ for all ${\rm W}\in(\mathbb{R}_{>0})^{2^{n}}$ and all $\sigma\in S_{n}$.

For example for any $1\leq i<j\leq n$, the permutation $(1i)(2j)\in S_{n}$ induces a bijection between the regions $O_{\{12\}}$ and $O_{\{ij\}}$ and hence

[TABLE]

Therefore, it suffices to consider one representative per $S_{n}$-orbit of regions. However, the number of required regions remains infeasibly large: We can identify each region $\mathcal{O}_{\mathcal{I}}$ with a hypergraph with vertices $[n]$ and edges $\mathcal{I}$ (a hypergraph is a generalization of a graph in which an edge can join any number of vertices). The number of required regions then coincides with the number of such hypergraphs up to isomorphism. Included in the hypergraphs of the regions are all graphs with vertices $[n]$, and the number of graphs on $n$ vertices up to isomorphism is well known to grow faster than any exponential function [GYZ14, §6.4].

Due to this problem, we preprocess the optimization problem. Note that we can rewrite Problem (1) in the following minimax form, as on each region $\mathcal{O}_{\mathcal{I}}$ the function $f_{\mathcal{I}}(W)$, as defined in (2), dominates all other $f_{\mathcal{I}^{\prime}}(W)$ with $\mathcal{I}^{\prime}\neq\mathcal{I}$:

[TABLE]

This can then be lifted to a problem with linear objective function and non-linear constraints:

[TABLE]

Given a fixed binding polynomial $P$, the polynomial system $\Phi(W)=P$ poses a serious computational problem. To sidestep this issue we now introduce new coordinates $s_{I}:=\prod_{I^{\prime}\subseteq I}w_{I^{\prime}}$. Here, $s_{I}$ represents the free energy difference of microstate $I$ to the fully unoccupied state. On the positive orthant $(\mathbb{R}_{>0})^{2^{n}}$, this defines a bijection

(\mathbb{R}_{>0})^{2^{n}}$$(\mathbb{R}_{>0})^{2^{n}}$$(w_{I})$$W=$$\Big{(}\displaystyle\prod_{I^{\prime}\subseteq I}w_{I^{\prime}}\Big{)}$$(s_{I})$$=S$$\Big{(}\displaystyle\prod_{I^{\prime}\subseteq I}s_{I^{\prime}}^{(-1)^{|I\setminus I^{\prime}|}}\Big{)}$$\varphi$$\phantom{{}^{-1}}\varphi^{-1}$$\in$$\in$$\in$$\in

In the new coordinates, the formerly polynomial constraints $a_{k}=\sum_{|I|=k}\prod_{I^{\prime}\subseteq I}w_{I^{\prime}}$ are simplified to linear constraints $a_{k}=\sum_{|I|=k}s_{I}$ for $k=1,\ldots,n$. However, the functions $f_{\mathcal{I}}$ become more complicated. For example, for $n=4$ and $\mathcal{I}=\{\{1,2,3\}\}$

[TABLE]

Since the $f_{\mathcal{I}}\circ\varphi^{-1}(S)$ remain monomials (with possibly negative exponents), this complication is of little consequence for the resulting problem:

[TABLE]

Example 3.2

Consider the polynomial $P_{1}$ of Example 2.7. We use SCIP [Gle+18], which is currently one of the fastest solvers for non-linear programming, to solve the resulting polynomial optimization problem (5). It uses a branch and bound method to solve the optimization problem with non-linear constraints.

Figure 3 shows the full input on the left and the partial output on the right. In the input, the first constraints c1, c2, c3 enforce $\sum_{|I|=k}s_{I}=a_{k}$ for $i=1,2,3$. The remaining constraints c4 to c19 enforce $f_{\mathcal{I}}\circ\varphi^{-1}(S)\leq r$. For example, c19 states that $s_{1}s_{2}s_{3}/s_{123}\leq r$ which is equivalent to $f_{\emptyset}(W)=(w_{123}w_{12}w_{13}w_{23})^{-1}\leq r$ in the coordinates $w_{I}$. The output states an approximate optimal value of $y=13.5$ and clearly lists all values of $s_{I}$ of the minimal molecules, which gave rise to the values of $w_{I}$ in Example 2.7.

Since the number of constraints are not too many when $n=3$, SCIP works well. However, for $n=4$ there are $2^{11}$ constraints which makes it a hard problem to solve directly.

To simplify our problem computationally, we compute an upper bound $b_{+}$ and a lower bound $b_{-}$ for the minimal absolute interaction of a given binding polynomial. If both bounds are identical, the minimum is found. To compute the upper bound we minimize the minimal absolute interaction in a single region. For the experiments in Subsection 4 we consider the orthant $\mathcal{O}_{\emptyset}$. The result may not be the global minimum - since the minimum may lie in another orthant - but it gives an upper bound. For the lower bound, we consider the relaxed problem by looking at a problem with less constraints. We pick a representative for each of the orbit of the group action $S^{n}$ and consider the corresponding region. We then only consider the constraints given by functions $f_{\mathcal{I}}$ dominant in those regions. Considering less constraints means that the feasible set of the optimization problem is bigger than originally allowed. Thus, this new problem gives the lower bound on the true optimization problem. When the gap between upper and lower bound is sufficiently small, we assume to have found a molecule minimizing the problem.

4. Experimental Results

In this section, we run computational experiments on data from [Ima73, Ike+83] which is also summarized in [CRBG86].

4.1. Data

The first data set consists of eight oxygen binding polynomials of human adult hemoglobin in its native form and of three chemically treated versions. The oxygen binding to each of the four hemoglobin variants has been observed under two different environmental conditions [Ima73]. The coefficients $a_{1}$, $a_{2}$ and $a_{3}$ of the corresponding polynomials $P_{1}$ to $P_{8}$ are listed in Figure 6. The coefficients $a_{0}$ and $a_{4}$ are standardized to $1$. Figure 4 describes the experimental conditions under which binding polynomials $P_{1}$ to $P_{8}$ were obtained and illustrates the ranking of the degree of cooperativity as defined by the maximal slope of the Hill plot $n_{max}$ (as reported by [Ima73]). The slope to the Hill plot relates the variance of the probability distribution on the macrostates $\{0,1,2,3,4\}$ at the respective ligand activity to the variance of a binomial distribution with the same mean [Abe05, MDH16]. A value larger than 1 indicates a variance higher than the maximal value an independent system can generate.

The second data set [Ike+83] contains five binding polynomials ($P_{9}$ to $P_{13}$) of native hemoglobin HbII of Scapharca inaequivalvis measured at different temperatures (Figure 5) [Ike+83, Table 3].

4.2. Minimal absolute interactions

Unlike the case $n=3$ in Example 2.7, solving Problem (5) for $n=4$ proved to be too hard for SCIP, which is why we list upper ($b_{+}$) and lower ($b_{-}$) bounds for the minimal absolute interaction $\|P_{i}\|$ instead. The upper bound is the minimal absolute interaction inside the region $\mathcal{O}_{\emptyset}$, i.e., it is obtained by solving Problem (3) for $\mathcal{I}=\emptyset$. In other words, the upper bound represents the minimal value on the region in which all interactions are equal to or larger than 1. The lower bound is computed through a relaxation of Problem (5). Inspired by the action of the symmetric group $S_{n}$ on the regions $\mathcal{O}_{\mathcal{I}}$, we picked a set of $180$ representatives and discarded all constraints $f_{\mathcal{I}}(W)\leq r$ with $\mathcal{I}$ not in them, see Section 6.

Figure 6 illustrates that for most binding polynomials our upper and lower bound allows for a sufficient approximation of the absolute minimal interaction. However, there remains a nontrivial gap between both bounds for $P_{2}$, $P_{3}$, $P_{4}$. Let us assume that the correct value for $P_{2}$ is close to the upper bound which we determined. Then, nearly all relations of the degree of cooperativity of polynomials $P_{1}$ to $P_{8}$ described by the maximal slope of the Hill plot $n_{max}$ (Figure 4) are also found with the minimal interaction criterion. The only exception is the ranking of $P_{4}$ and $P_{6}$. Whereas $P_{4}$ describes a more cooperative system than $P_{6}$ according to $n_{max}$, we obtain a different picture when using the minimal interaction criterion which states that $P_{6}$ is almost as cooperative as $P_{2}$. The question arises why the two measures of cooperativity are so different. An answer may be found considering the coefficients of $P_{4}$ and $P_{6}$. The coefficient $a_{2}$ of $P_{6}$ is larger than $a_{1}$ and $a_{3}$, which indicates that the macrostate of having $2$ sites occupied –at sufficiently high ligand activity– is more probable than the states of having $1$ or $3$ ligands bound. To generate this type of distribution, we may need absolute interaction energy at each macrosotate level, for instance negative and positive interaction energy. Moreover, a family of distributions with a relatively high weight on macrostate $2$ has a reduced possibility to generate a high variance since an extraordinary high variance requires more weights on the macrostates [math] and $4$.

The polynomials $P_{9}$ to $P_{13}$ were derived from the same protein at different temperatures [Ike+83]. We see that the reported $n_{max}$ is not behaving monotonously with increasing temperature, but that it varies between $2.08$ and $2.12$. Since one might expect a more monotonous development, we interpret this observation as a result of measuring imprecision of the coefficients. Comparing the minimal interaction to $n_{max}$, we see the same problem on a larger scale, which is a results of the minimal energy criterion being defined on an exponential scale of actual physical energies.

We consider the molecules realizing the upper bound in more detail.

4.3. Molecules solving the optimization problem

Figure 7 presents the molecules realizing the upper bound of the minimal absolute interaction for all polynomials. We see that the major part of the estimated interactions is very close to $1$. This is not surprising since the considered optimization problem minimizes the absolute values of the physical energies required to generate the binding curve. $L_{1}$ optimizations –for instance such as the least absolute shrinkage and selection operator (Lasso) regression [Tib96]– are known to provide sparse solutions, meaning that many of the physical energies of the solution will be zero. Thus, this observation fits well to the context. Considering the structure of the solutions in more detail, we see that most solutions include a non-trivial weight for a second degree interaction ($w_{34}$). The only polynomial whose minimizing molecule in $\mathcal{O}_{\emptyset}$ does not possess an interaction of pair-wise interaction is $P_{2}$ which has the highest $n_{max}$. This observation seems plausible since the weights should be more towards the extreme values, that is to the macrostates [math] and $4$ to obtain a high variance of the macroscopic distribution. Compared to $P_{1}$, $P_{2}$ puts more probability on the macrostates $3$ and $4$. With this view, we can also explain the apparent contradiction that we find when using $n_{max}$ or the minimal absolute interaction to compare polynomials $P_{4}$ and $P_{6}$. As described above, $P_{6}$ has a coefficient $a_{2}$ which is larger than $a_{1}$ and $a_{3}$ and almost of same size as $a_{4}$, which is $1$. This means that the distribution on the macrostates is stabilized towards the macrostate $2$, which hinders very large variances and thus leads to a smaller $n_{max}$ than that of $P_{4}$. However, this circumstance leads to a high interaction weight for a microstate of macrostate $2$. While the main weight for $P_{4}$ is on macrostate $4$ -which helps to achieve a higher value of $n_{max}$, $P_{6}$ requires more interaction energy to model the behavior at macrostate $2$.

5. Discussion and Outlook

5.1. Cooperativity and minimal absolute interaction

A commonly used macroscopic conceptualization of cooperativity is that “It’s all or nothing” [HA09], meaning for hemoglobin that the probability weights of the extreme macrostates [math] and $4$ are relatively big. This definition also aligns well with using the maximal Hill slope as criterion of cooperativity, since it detects a variance which is higher than that of any independent system. However, observations made on the macroscopic level have often been interpreted with a view on microscopic properties, for instance that “binding of a ligand molecule increases the receptor’s apparent affinity, and hence increases the chance of another ligand molecule binding” [SL13]. Connecting macroscopic and microscopic properties is a non-trivial problem if no additional restriction on the underlying molecule are made. Microscopic properties can be easily inferred from the macroscopic behavior if additional assumptions are made, for instance the molecule being composed of physically indistinguishable sites. In particular, given a binding polynomial, we can calculate all interaction energies easily when we assume that $w_{ij}=w_{kl}$, $w_{ijk}=w_{lmn}$, etc. However, without additional knowledge or restrictions, it is difficult to relate macroscopic and microscopic properties of larger, asymmetric systems [MDH16]. The minimal interaction energy required to generate the observed binding curve is a good candidate to connect macroscopic and microscopic characteristics and additional side conditions based on knowledge of the molecule structure could also be incorporated [Mar17].

As illustrated by the examples of $P_{4}$ and $P_{6}$, the minimal interaction energy criterion captures a conceptually different type of cooperativity than $n_{max}$. The polynomial $P_{6}$ requires a relatively high interaction energy due to giving the highest weights to macrostates [math], $2$ and $4$, but stabilizing towards the intermediate value $2$ prevents an extreme variance. Contrary, $P_{4}$ has monotonously decreasing coefficients $a_{1}$ to $a_{3}$. Here, modeling the distribution requires less interaction energy to capture the relation of the probability on macrostates $1$, $2$ and $3$, and a higher value of the variance of the distribution can be realized.

5.2. The mathematical problem

Ligand binding raises many intriguing questions in applied algebra and geometry, and stands to benefit equally from the provided insight into the intrinsic geometry of the problem. While previous works have explored the application of polynomial system solving techniques to a situation with different types of ligands [RMT19, MSU13], this work explores the use of polynomial optimization. Using the minimal interaction energy as a criterion for cooperativity without assuming the target molecule to be composed of physically identical units, we deal with an algebraic optimization problem. Even though there exists dedicated software to solve this kind of problem, we had to apply a transformation of the coordiantes (from $w_{I}$ to $s_{I}$), a transformation of the optimization problem (Eq.(4)) and approximations by upper and lower bounds to approach the considered examples. This effort illustrates that the setup of looking for the minimal absolute interaction by this method becomes rapidly complicated with increasing number of binding sites. Whereas the case of three binding sites can be solved easily with SCIP, tetrameric hemoglobin with four binding sites already posed computational problems. A way to facilitate the problem may be to exploit its intrinsic structure. For instance, there are results for strongly symmetric systems which might be generalized and tailored to the system with which we are dealing [RS18]. In addition, other mathematical fields such as tropical geometry [MS15] may provide interesting tools to approach this problem.

5.3. Other mathematical perspectives

5.3.1. Tensors

Tensors [Lan12] are multidimensional arrays of real numbers which generalize matrices. We can regard molecules as (positive) tensors in $\mathbb{R}^{2\times\dots\times 2}=(\mathbb{R}^{2})^{\otimes n}$, whose entries are indexed by $\{0,1\}^{n}$ which is in bijection with the subsets of $[n]$.

From this point of view, many concepts in tensor theory have direct analogues in ligand binding, such as symmetric tensors and molecules with indistinguishable sites. Notably, the central concepts of tensor rank and tensor rank decomposition have a special meaning in ligand binding:

Proposition 5.1

Let $S=(s_{I})_{I\subseteq[n]}\in\mathbb{R}^{2^{n}}$ be a molecule with $n$ sites and whose corresponding tensor $(s_{I})_{I\in\{0,1\}^{n}}$ has (positive) rank $k$. Then $S$ can be written as a weighted average of $k$ molecules with $n$ independent sites. In particular, a corresponding mixture of the $k$ molecules will have the same binding curve as $S$. Moreover, this is not possible with $k-1$ molecules with $n$ independent sites.

Proof.

Rank-one tensors are of the form $v_{1}\otimes\dots\otimes v_{n}$ for some $v_{i}=(s_{i},t_{i})\in\mathbb{R}^{2}$ with entries

[TABLE]

They can be regarded as molecules with particularly simple interaction, as molecules with independent sites are rank-one tensors with $t_{i}=1$. A tensor $S$ of rank $k$ can be expressed as

[TABLE]

for $k$ rank-one tensors $v_{1}^{(j)}\otimes\dots\otimes v_{n}^{(j)}$ and cannot be expressed as a sum of $k-1$ rank-one tensors. Supposing $v_{i}^{(j)}=(s_{i}^{(j)},t_{i}^{(j)})\in\mathbb{R}^{2}$, we can reformulate this into

[TABLE]

thus $S$ is a weighted average $k$ molecules with $n$ independent sites. ∎

5.3.2. Tropical Geometry

Tropical geometry [MS15] revolves around functions on the max-plus semiring. One example of such as function is the absolute interaction in logarithmic coordinates:

(\mathbb{R}_{>0})^{2^{n}}$$\mathbb{R}^{2^{n}}$$(\log(w_{I}))$$=:(x_{I})$$(w_{I})$$\displaystyle\prod_{I\subseteq[n]}\max(w_{I},w_{I}^{-1})$$\displaystyle\sum_{I\subseteq[n]}\max(x_{I},-x_{I})$$=:\|(x_{I})\|$$\in$$\in

In tropical geometry, the function $\|\cdot\|$ is referred to as a tropical rational function. It is a piecewise-linear function, linear on the image of each $\mathcal{O}_{\mathcal{I}}$, and the boundary separating its regions of linearity is called its tropical hypersurface $\text{Trop}(\|\cdot\|)$. The gradient of $\|\cdot\|$ is constant on each region of linearity, which is why understanding the intersection of a set $L\subseteq\mathbb{R}^{2^{n}}$ with $\text{Trop}(\|\cdot\|)$ is paramount to computing the minimum of $\|\cdot\|$ on $L$.

In our setting, $L$ is the image of the feasible set given by the conditions $\Phi(W)=P$ which is an exponential variety of codimension $n$. Since the number of regions of linearity is finite, there exist only a finite number of intersection patterns of $L$ and $\text{Trop}(\|\cdot\|)$. Classifying these intersection patterns depending on the coefficients of $P$ will contribute greatly towards computing the minimal interactions and minimizing molecules.

5.3.3. Dependency measures

From a stochastic or statistical point of view, the problem of cooperativity is a problem of quantifying (minimal) dependency, which is a classical research topic in these fields [Rén59, SW81, Koy+87]. In particular copulas [Sch91, DS10, Nel07] have become a common tool to describe dependencies between random variables. Copulas cannot be applied directly to the problem since they are defined on continuous, not discrete, random variables. Moreover, they aim at modeling the dependency of given random variables which in our scenario corresponds to the dependency of the sites of a given molecule. The minimization problem of finding the minimal dependency required to generate a given distribution of a sum of dependent variables is usually not included. Nevertheless, the problem of cooperativity should also be treated from a more probabilistic point of view and the available concepts should be considered in more detail in the context of ligand binding.

5.4. Conclusion

We illustrated how searching for the minimal interaction energy required to generate a binding curve leads to a problem of algebraic optimization. Even though, there are computational tools available and the structure of the problem seems relatively simple on first sight, it becomes complicated when more than $3$ binding sites are considered. In particular the fact that we had to calculate a lower bound by a relaxed problem and not having been able to close the gap between upper and lower bound completely, indicates that our approach can currently not be enveloped in a simple ready-to-use tool to quantify cooperativity. We showed that the concept of minimally present interaction energy is different from using the maximal Hill slope as a measure of cooperativity. The first also detects required changes in affinity from one macrostate to the other, whereas the latter only focuses on the macroscopic variance of the distribution. Future work may use other mathematical concepts such as tensor theory, tropical geometry and probabilistic measures of dependencies to approach questions related to cooperativity.

6. List of Constraints

Below is a compact representation of the list of $180$ constraints we considered to compute the lower bounds $b_{-}$ of minimal absolute interaction in Figure 6. There, $\{12,34,123\}$ represents the constraint $f_{\mathcal{I}}\leq r$ for $\mathcal{I}=\{12,34,123\}$, which is the absolute interaction $\|W\|$ for $W\in\mathcal{O}_{\mathcal{I}}$ where $w_{12},w_{34},w_{123}\leq 1$ and $w_{I}\geq 1$ otherwise.

$\emptyset$, $\!\{12\}$, $\!\{123\}$, $\!\{1234\}$, $\!\{12$, $\!\!13\}$, $\!\{12$, $\!\!34\}$, $\!\{12$, $\!\!123\}$, $\!\{12$, $\!\!134\}$, $\!\{123$, $\!\!124\}$, $\!\{12$, $\!\!1234\}$, $\!\{123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14\}$, $\!\{12$, $\!\!13$, $\!\!23\}$, $\!\{12$, $\!\!13$, $\!\!24\}$, $\!\{12$, $\!\!13$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!124\}$, $\!\!\!\{12$, $\!\!34$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!234\}$, $\!\{12$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!123$, $\!\!134\}$, $\!\{12$, $\!\!134$, $\!\!234\}$, $\!\{123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!1234\}$, $\!\{12$, $\!\!34$, $\!\!1234\}$, $\!\{12$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!134$, $\!\!1234\}$, $\!\{123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!134\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!124$, $\!\!234\}$, $\!\{12$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!123$, $\!\!134$, $\!\!234\}$, $\!\{123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\!\!\{12$, $\!\!123$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124$, $\!\!234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!124$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!123$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$, $\!\{12$, $\!\!13$, $\!\!14$, $\!\!23$, $\!\!24$, $\!\!34$, $\!\!123$, $\!\!124$, $\!\!134$, $\!\!234$, $\!\!1234\}$

Acknowledgments. Authors thank Felipe Serrano for valuable technical support with SCIP